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ABSTRACT 


MEASUREMENT SYSTEM CHARACTERIZATION IN THE 
PRESENCE OF MEASUREMENT ERRORS 


Sean A. Commo 
Old Dominion University, 2012 
Director: Dr. Drew Landman 


In the calibration of a measurement system, data are collected in order to esti- 
mate a mathematical model between one or more factors of interest and a response. 
Ordinary least squares is a method employed to estimate the regression coefficients 
in the model. The method assumes that the factors are known without error; yet, 
it is implicitly known that the factors contain some uncertainty. In the literature, 
this uncertainty is known as measurement error. The measurement error affects both 
the estimates of the model coefficients and the prediction, or residual, errors. There 
are some methods, such as orthogonal least squares, that are employed in situations 
where measurement errors exist, but these methods do not directly incorporate the 
magnitude of the measurement errors. This research proposes a new method, known 
as modified least squares, that combines the principles of least squares with knowl- 
edge about the measurement errors. This knowledge is expressed in terms of the 
variance ratio - the ratio of response error variance to measurement error variance. 
The variance ratio takes on values between 0 and 1, and for calibration applications, 
the ratio is typically less than 0.0625. In addition to modified least squares, a new 
definition of residual errors based on the variance ratio is proposed. Through several 
simulation studies, it is observed that the new estimator can yield different esti- 
mates of the regression coefficients and improve the residual error over ordinary least 
squares. As a result, modified least squares is shown to be an alternative estimation 
method in the presence of measurement errors. 
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NOMENCLATURE 

a = Distance Minimized by Ordinary Least Squares 
d = Distance Minimized by Modified and Orthogonal Least Squares 
E = Expected Value Operator 
EFT-1 = Exploration Flight Test-1 
F.S. = Full-Scale 

k = Number of Factors 
ME = Measurement Error 
MLS = Modified Least Squares 
MSE = Mean Squared Error 
n = Number of Design Points 
N = Normally Distributed 

NIST = National Institute of Standards and Technology 
OLS = Ordinary Least Squares 
OrthLS = Orthogonal Least Squares 

p = Number of Model Terms without the Intercept 
S xx = Sum of Squares in x 
S X y = Sum of Cross-Products in x and y 
S yy = Sum of Squares in y 
u = Random Error or Measurement Error in x 
Var = Variance Operator 


W = Factor with Error 
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x = Factor without Error 
x = Sample Mean of x 
x = Model Matrix 

y = Observed Response without Error 

Y = Observed Response with Error 

Y = Sample Mean of Y 

Y = Predicted Response of Y 

a = Angle Formed by d and Estimated Model 
/ 3 = Regression Coefficient 
j3 = Estimated Regression Coefficient 
e = Random Error in y 
7 = Variance Ratio 

c f> = Angle Formed by a and Estimated Model 
a = Estimated Measurement System Accuracy 
cr 2 a = Variance of u 
cr e = Standard Deviation of e 
a \ = Variance of e 

E = Variance-Covariance Matrix of the Measurement Errors 
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CHAPTER 1 


INTRODUCTION 


1.1 BACKGROUND 

Regression analysis is a collection of statistical methods and tools used to es- 
timate mathematical relationships between one or more explanatory variables, or 
factors, and responses. Regression can be classified into several categories includ- 
ing simple or multiple, linear or nonlinear, and parametric or nonparametric. The 
earliest published work on regression analysis was done by Legendre and Gauss in 
the early 19th century although Gauss began formulating his ideas in the late 18th 
century (Draper and Smith, 1998). Both Legendre and Gauss independently derived 
what is now known as the method of least squares and applied their methods to 
orbital mechanics. Additionally, Gauss developed the Gauss-Markov theorem, which 
revealed an elegant property of the least squares method. The method of least squares 
minimizes the sum of the squares of the errors. The least squares estimator is the 
best linear unbiased estimator (BLUE). In other words, the regression coefficients 
estimated from least squares have minimal variance and are unbiased. While often 
overlooked, this result is useful when applying the least squares method. 

As Legendre and Gauss first recognized, the rationale behind estimating these 
mathematical relationships ranges from validating physical laws and phenomenon to 
understanding system behavior and performance. The relationship representing a 
physical law or system response in terms of k factors is expressed as 


y = f(x x k ). 
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In fc-dirnensional space, the function f (x i, . . . ,Xk ) is known as a response surface. 
The functions that make up the true response surface may be simple or complex. 
Because the true functions are seldom known, a simpler function over a small region 
of interest is used to approximate f (x i, . . . ,Xk)- The approximate relationship can 
be expressed as 

y = g(xi,-..,x k ) + e (l) 

where e is the error. For this research, the form of g (aq, . . . ,Xk ) is limited to a class 
of linear models based on a Taylor-series expansion in k factors. These functions are 
often sufficient in approximating the true function since they are extremely flexible. 
Furthermore, estimating the coefficients in the function from historical or experi- 
mental data is accomplished using the method of ordinary least squares (OLS). In 
the latter case, experiments are designed to collect sufficient information in order to 
estimate the assumed form of g (aq, . . . , aq). 

Second-order Taylor-series response surface functions often work well in real- world 
applications and therefore are one of the most commonly estimated functions (Myers 
et al., 2009). The mathematical form of this model is 

k k— 1 k k 

y = (3o + ^2 @ iXi + SI PiJ x i x i + + e ( 2 ) 

2=1 2—1 j=i-\-l 2=1 

where k is the number of factors, the /3’s are the regression coefficients, and e is the 
error. Based on Equation (2), there are 


1 + 2 k + 


k(k- 1) 
2 


1 +p 


regression coefficients in the model. As a requirement, the experimental design must 
contain 1 + p unique design points. If x is the model matrix which contains the 
design points, then the experimental design is said to have a sufficient number of 
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unique design points when (x 1 x) is full rank. Additionally, because the function is 
second-order, the experimental design must also contain 3 unique levels of each of 
the k factors. 

The use of OLS to estimate the regression coefficients carries the following as- 
sumptions about the factors, responses, and mathematical model: 

• Appropriate model specification. The form of the mathematical model is lin- 
ear with respect to the regression coefficients, and the model lack-of-fit is not 
significant. 

• Linear independence. The factors in the model matrix, x, are linearly inde- 
pendent or the responses are a linear combination of the factors and regression 
coefficients. From the experimental design perspective, linear independence is 
obtained through a sufficient number of unique design points to fit the assumed 
mathematical model. Linear independence is expressed mathematically when 
x is of full column rank or (x 3 x) is well-conditioned. 

• Independent errors. The errors, e, across all the responses are independent. 
Mathematically, this is expressed as 


Cov (e;, 6j) = 0 for i ^ j. 

For errors that are not independent, generalized least squares (GLS) should be 
used. 

• Homoscedasticity or constant variance. The variance of the errors, Var(e), 
across all the responses are equal. This is independent of the settings of the 
factors. 

• Weak exogeneity. The factors are fixed variables and therefore assumed to 
be known without error. From the regression model, the errors across all the 
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responses have a conditional mean of zero and are independent of the factors. 

Other assumptions about normality and identically independent errors are not nec- 
essary for OLS, but these assumptions do support additional properties of the esti- 
mators. 

The error, e, in the response surface function is attributed to one or more of the 
following: 

• The response, y, is observed with error, which may result from systematic biases 
or random fluctuations. This is the most commonly assumed cause of error in 
the estimated model. 

• The form of the response surface function, g (xi, . . . ,Xk), is incorrect due to: 

— factors that were excluded but affect the response, or 

— a more complex true function, f (x i, . . . ,Xk), than the assumed function. 

• The factors are not known without error. Typically, any error in the factors 
is neglected so the OLS can be used to estimate the mathematical model. 
In the literature, errors in the factors are referred to as errors-in-variables or 
measurement error (ME). This is the motivation behind the research presented 
in this dissertation. 

The source of MEs can be attributed to one or more causes, including instrument 
and sampling errors. In practice, all physical experiments contain some degree of 
ME, but MEs are usually considered negligible in favor of using OLS to estimate a 
model. Depending on the application, the MEs may be small relative to the other 
sources of error in the mathematical models and therefore are assumed not to impact 
the modeling. However, there are a few recommendations for how small the MEs 
should be in order to neglect them in the regression analysis. Furthermore, statistical 
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software packages that recognize MEs in the modeling are extremely limited, typically 
limited to the simple linear model or 


y = A) + Pix + e. 


At the National Aeronautics and Space Administration (NASA), MEs occur in many 
complex systems and applications, and the simple linear model is not adequate. For 
example, consider a wind-tunnel experiment where the factors of interest are angle- 
of-attack, yaw angle, and Mach number and the responses are aerodynamic forces 
and moments, and the objective of the experiment is to estimate the mathematical 
relationships between the factors and responses. It is known implicitly that both the 
factors and the responses contain some uncertainty; yet, there is no formal method- 
ology employed to estimate the response surface function between the factors and 
responses which considers both sources of error. This issue is further exemplified 
when the same test article is tested in different facilities and differences in the mod- 
eled relationships cannot be reconciled. Since the ultimate goal is an understanding 
of the true relationship between the factors and responses, any unaccounted for errors 
can complicate the ability to achieve this objective. 

The research presented in this dissertation proposes a new general methodology 
for mathematical modeling in the presence of MEs. This methodology is applicable 
to simple and complex systems and enables direct incorporation of uncertainties in 
all variables, which is often not considered in practice. 

1.2 APPLICATIONS 

Both within and outside of NASA, there exist several examples of ME. At NASA, 
complex measurement systems are designed to collect data for various projects across 
all the Mission Directorates, including Aeronautics, Exploration, and Atmospheric 
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Science. In complex, custom-designed measurement systems, system-level calibra- 
tions are critical in understanding and capturing the true performance. However, 
there exists a significant difference between classical instrument calibration and the 
characterization of a complex measurement system. In laboratory calibration, the 
factors are considered to be known without error and are traceable to standards 
from the National Institute of Standards and Technology (NIST) within the United 
States. In general, many of NASA’s measurement systems cannot be calibrated with 
a system-level standard that is traceable to NIST. 

In aeronautics, scaled testing of flight vehicles is performed in wind tunnels to un- 
derstand aerodynamic performance. The primary measurement system used during 
these wind-tunnel experiments is a force-balance. A force-balance is a multiple-axis 
load cell that provides simultaneous, high-precision measurements of aerodynamic 
forces and moments in up to 6 degrees of freedom. The number of aerodynamic 
components that can be measured by a force-balance is a function of its mechanical 
and electrical design. Through a calibration experiment, a mathematical model can 
be developed between the applied forces and moments, or factors, and the electrical 
outputs of the strain gauges, or responses. The literature provides several examples 
of improvements in calibration techniques, including calibration in the presence of 
additional factors, such as temperature (Lynn et al., 2012; Parker et al., 2001; Parker 
and DeLoach, 2001). However, current force-balance calibration methods rely on the 
assumption that the applied forces and moments are known without error. For ex- 
ample, in force-balance calibrations that use gravity-based loads, or dead-weights, 
the individual weights are calibrated against NIST standards; therefore, each weight 
is known within a small uncertainty. These uncertainties propagate into errors in the 
applied forces and moments. 

In exploration, research in the development of newer, more robust entry, descent, 
and landing (EDL) technologies is a major focus. Both landed payload mass and 
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landing accuracy are driven by the EDL architecture for a given mission. The cur- 
rent unmanned exploration mission to Mars, the Mars Science Laboratory (MSL), 
is nearing the limit of the capabilities of available EDL technologies. While un- 
manned missions have less stringent requirements for landing, manned missions are 
infeasible with the current EDL systems (Braun and Manning, 2006). The MSL 
Entry, Descent, and Landing Instrumentation (MEDLI) was proposed to address 
some of the challenges associated with the development of newer, more robust EDL 
technologies (Gazarik et al., 2008). MEDLI is a suite of sensors installed on the 
forebody heatshield of the MSL entry vehicle. One of the subsystems of MEDLI 
is the pressure measurement system, which consists of two components: the Mars 
Entry Atmospheric Data System (MEADS) and the associated electronics. MEADS 
is a series of through-holes, or pressure ports, in the heatshield that connect via 
stainless steel tubing to pressure transducers. Power is provided to MEADS by the 
electronics system. Commo and Parker (2012) discuss a system-level calibration ap- 
proach that was developed and employed to characterize the performance of the pres- 
sure measurement system. Typically for spaceflight instrumentation, component-by- 
component calibrations are performed with some of the components being calibrated 
with NIST-traceable standards. However, a system-level, NIST-traceable standard 
for the pressure measurement system does not exist. 

In atmospheric sciences, Earth’s upper atmosphere is receiving increased research 
emphasis. Researchers use satellite observations to gain insight into trends and causes 
of changes in the Earth’s climate. Missions such as the Clouds and Earth’s Radiant 
Energy System (CERES) have been proposed to answer questions about how radia- 
tive energy in the atmosphere affects climate change. This includes how increases in 
quantities of carbon dioxide affect the balance of energy. Satellite missions employ a 
series of simple instruments and complex measurement systems that are integrated 
to collect the appropriate data in order to support the mission objectives. These 



instruments and measurement systems are calibrated on the ground. However, satel- 
lites operate in a vacuum environment and are exposed to radiation; both of these 
environmental factors are known to influence the measurements. In-situ adjustments 
are made to the instruments and measurement systems in flight to correct for en- 
vironmental effects. Measurement standards are unavailable for on-orbit satellites, 
so any in-situ adjustments are relative to a source, which contains some error. Sim- 
ilar to the force-balance example, the uncertainty from the source propagates into 
measurements made by the sensors on the satellites. 

1.3 CLASSIFICATION OF MEASUREMENT ERRORS 

Since both fixed and random variables are discussed within this research, it is 
important to distinguish each mathematically. For consistency, Buonaccorsi’s (2010) 
convention is used to identify these variables. In the case of a fixed variable, a lower 
case letter (e.g. x) is used while a capital letter (e.g. X) indicates a random variable. 
Bold letters (e.g. x and X) are used to indicate a vector or matrix of fixed and random 
variables, respectively. The given operator is also widely used and is designated by 
“|". For example, \x can be interpreted as “given x.” When referring to the random 
variable X, then |x means “given X equal to x." Furthermore, distinguishing the 
difference between fixed and random factors is also important. If the factor x is taken 
to be fixed, then the model is identified as a functional model. The random factor 
X is the structural model. 

To demonstrate the difference between functional and structural models, consider 
an experiment where a factor is applied load. Suppose a 100 Ibf. load is to be applied 
and there are 30 possible configurations to apply the 100 Ibf. load. If a configuration 
is randomly selected for every instance that the load is applied, then the structural 
model is appropriate. However, if the same configuration is selected for every instance 
the load is applied, then the functional model is correct. For this application, the 



9 


focus is on the functional model, but many of the conclusions from the functional 
model can easily be extended to the structural model (Gleser, 1983). Both functional 
and structural models are classified as classical ME models, where W is the error- 
prone measurement of the factor x or X. 

1.4 GENERAL FORM OF THE MEASUREMENT ERROR MODEL 

The general form of the classical linear model is 


V — A) + (3iXi + . . . + f3 p Xp + e, 


where the (3s are the regression coefficients and e is the normally distributed random 
error of y with a mean of zero and a constant variance of a\. The subscript p is 
the number of model terms without the intercept and does not necessarily equal the 
number of factors in an experiment. For a second-order response surface model based 
on a Taylor series expansion, 


p = 2k + 


k(k- 1) 
2 


For example, if k — 2, then p = 5 and the model is 


y = (3 o + (3iXi + (3 2 x 2 + (3 12 x 4 x 2 + (3iix\ + (3 2 2 x\ + e 

= A) + (3ix i + (3 2 x 2 + (3 3 x 3 + (3 4 X4, + (3 5 x 5 + e, 


where f3 3 = /3 i2 , /ffi = /3 U , (3 5 = (3 2 2 , x 3 = xix 2 , x 4 = x'f, and x 5 = x\. Using 
summation notation, the second-order Taylor series model is 

k k— 1 k k 

y = P 0 + ^ PaXa + ftab%a%b ^ ^ Paa^a 

a = 1 a— 1 b=a+l a= 1 
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k k — 1 k 

where PaX a are the linear terms, PabX a Xb are the two-factor interaction 

a= 1 a= 1 b=a-\-l 

k 

terms, and Pa a x 2 a are the second-order terms. 

a — 1 

From an experiment, n observations are used to estimate the regression coeffi- 
cients in the model. At the ith observation, the classical model is 


Vi = Po + Pixit + • • • + P P x Pi + 6i i = 1, 2, . , . , n. 


For each observation, the expected mean and variance operators are applied to get 


E (yi \x u , . . . , x Pi ) = P 0 + PiXii + • . . + P P x Pi or E (e*) = 0, and 


Var {yi\x l0 x Pi ) = Var (e*) = of. 


As expressed, the variance is assumed to be constant across all observations. Fur- 
thermore, since the observations are assumed to be independently and identically 
distributed, the covariance between two responses, or Cov {yi,yj \x\ t , . . . ,x Pi ) where 
i ^ j, is zero. The n equations are expressed as 


y i — Po + Pix u + . . . + (3 p x pi + 6i 
V 2 = Po + P\X\ 2 + . . . + P p x P2 + 62 

2 /n = Po + PiX i„ + . . . + P p x Pn + e n 
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Combining the system of equations in matrix form yields 



y = x/3 + e 

where y is a n x 1 vector, x is a n x (p + 1) matrix, (3 is a (p + 1) x 1 vector, and e 
is a n x 1 vector. The expected mean and covariance, in matrix notation, are 

E (y) = E (x/3) + E (e) = x/3 or E (e) = 0, and 

Var (y) = Var (x/3) + Var (e) = of I. 

However, the factors are not known without error. Instead, the error-prone value of 
the factor x is 

Wi = Xi + Ui 

where iq is the random ME associated with x % and is distributed with a mean of 
zero and a variance of of . This form of ME is known as the additive model since 

lit 

E (Wi\xi) = Xi . Therefore, Mf is unbiased for the true value of x^. For k factors, the 
errors in the xs are 

((/],..., iif,-) ~ [ 0 , o 2 Uk )} . 

Constant variance of the errors in the xs is not necessary, but for most practical 
applications, this assumption of constant variance is sufficient. The xs in the classical 
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model are replaced with the ME model to get 


V — Po + Pi (^1 + u i ) + • • • + P p ( x p + Up) + e 


— (3q + + • • • + (3pW p + e. 


Unlike the classical model, the statistical properties of the ME model cannot be 
generalized. Starting with the first-order model (p = k ) 

k 

V = Po + PaWa + 6 , 

a=l 

the expected value and variance are (Buonaccorsi, 2010) 

k 

E (y\Wi, W k ) = p 0 + ^PaXa, and 

a— 1 
k 

VarW Wu---,W h ) = Y.0K. + <%- 

a= 1 

If the MEs are correlated, then the expected value and variance include the co- 
variances between u ±, ... ,u k . However, independence is assumed between the MEs. 
Regardless of independence, the variance of the response is inflated due to the ME. 
For a first-order with interaction model (p — k + [k (k — 1)] /2), 

k k — 1 k 

V — Aj + PaWa + EE P ab W a W b + e, 

a= 1 a= 1 b=a+ 1 

the expected value and variance are (Buonaccorsi, 2010) 

k k— 1 k 

E(y\W u ..., W t ) “A + E^ + E E PabX a x b , and 

a= 1 a=l fe=a+l 
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Vai{y\W 1 ,,..,W k ) = '52Pfc 


d=l 


k - 1 

+e 

a= 1 b=a+l 


E K( 


2 2 


, 2 2 
+ 


L ^2 2 


)] + 


O', 


Similar to the first-order model, the first-order with interaction model is additive 
since response is unbiased for independent MEs. The variance is slightly larger 
for the first-order with interaction model since it contains additional model terms. 
Lastly, for the second-order model (p = 2k + [k (. k — 1)] /2) 


k k — 1 k k 

y — A) + fiaWa + EE PabW a W b + J^PaaW^ + 6, 

a= 1 a= 1 b=a+ 1 a=l 


the expected value and variance are (Buonaccorsi, 2010) 

k k—1 k k 

E (y\w u .... W k ) = ft, + E /W + E E PabXaXb + Aia (a^ + , and 

a= 1 a = 1 b=a + 1 a=l 


Var ( 2 /| Wi, . . . , hFfc) 


a=l a=l 6=a+l 

+ X] [ 4 ^aa a: o 0 'u a +(T 1] + ' 

a=l 


Although the ME is additive for any single factor, the second-order model is nonad- 
ditive since E (y\Wi , . . . , 14 A) contains some bias. More specifically, the bias in the 

k 

second-order model is ^ (3aa&u a ■ The variance is further inflated by the second-order 

a= 1 

terms in the model. Therefore, the effect of ME on the statistical properties of the 
model increases as the order of the model increases. A summary of the bias and 
variance results are presented in Table 1. 



TABLE 1: Summary of Bias and Variance for Classical and Measurement Error Models 





15 


1.4.1 EXAMPLE: FORCE TRANSDUCER CALIBRATION 

Suppose the output of a force transducer, y, is a function of an applied force, 
F, and temperature, T. The applied forces and temperatures have been calibrated 
against known standards with the following uncertainties: 

uf ~ N (() lbs., 0.01 lbs. 2 ) , and 
ut ~ iV ( OK, 0.01 K 2 ) . 

In physical terms, temperature is known with a precision of VO.OI K 2 = 0.1 degrees 
Kelvin. If the true model for the force transducer in the presence of ME is 

V — A) + A. (F + uf) + /?2 (T + ut) + [3i2 (F + uf) (T + ut) + e 

where /3 0 = 0, = /?i 2 = 1, and e ~ N (0, 0.01), then the expected mean is 

E (y\F, T, uf, ut) = f3o + PiF + (3 2 T + /3 i2 (F x T) 

= F + T + (FxT). 


Additionally, the variance is 

+ '| + r2 + Pi2 (n n UT + T G Up + a Up a UT 

0.01 + 0.01 + 0.01 + (0.01T 2 + 0.01F 2 + 0.0001) 

0.0301 + 0.01T 2 + 0.01F 2 . 

In the classical linear model, the variance of y is 0.01, which comes from e. However, 
since variance is never less than zero, the variance in the classical model will always 
be an underestimate of the true, unknown variance. 


Var (y\F, T, u F , u T ) = 
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CHAPTER 2 


LITERATURE REVIEW 


Adcock (1877, 1878) is credited with developing the first least squares-based method 
that is applicable to measurement error (ME) models. He proposed a multi- 
dimensional objective function as an alternative to the one- dimensional objective 
function that is used in ordinary least squares (OLS) estimation. Adcock considered 
the normal distance from the estimated line to an observation. Since he utilized least 
squares principles to estimate the regression coefficients, Adcock minimized the sum 
of the squared normal distances over n observations. This approach is the basis of 
orthogonal least squares (OrthLS) regression. It is noted that the line estimated by 
OrthLS always passes through the centroid of the data. However, OrthLS is limited 
since it weights equally all the variables in the model. In terms of ME models, it is 
equivalent to assuming the variability in the factors is equal to the variability in the 
response. Meanwhile, Pearson (1901) showed that the line estimated by OrthLS is 
bound by the OLS regressions of Y on x and x on Y in two-dimensional space. 

Kummel (1879) expanded upon the work of Adcock by developing a more general 
solution that is valid when the errors are not assumed to be equal. He is the first 
person to define the concept of variance ratio, 7 ; mathematically, the variance ratio 
is expressed as 

7 = 4 - 

where of is the variance of the error in the response and of is the variance of the 
ME. The variance ratio is an important parameter later in this research as it is 
used to derive a new estimator. Based on prior experimental data and subject- 
matter expertise, Kummel argues that obtaining an estimate of 7 is reasonable in 
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practice. In most calibration equations, 7 takes on values between 0 and 0.0625. The 
interpretation of these values is discussed in the next chapter. 

Deming (1931, 1964) continued developing the work of Adcock and Kummcl and 
applied it to the special case of the simple linear model. The method he derived 
is appropriately named Deming regression. The error-prone observations in a set of 
experimental data are 

Yi = Vi + 7 


and 


Wi = Xi + Ui 


where i — 1, 2, . . . , n. The errors e and u are independent random variables and are 
assumed to be distributed with a mean of zero and a constant variance of of and of, 
respectively, and are related through the variance ratio, 7 . If the variances are not 
constant, the method is invalid. The model coefficients are estimated in Dcming’s 
method by minimizing 

n / 2 2 \ n 

E \% + % ) = E [« - ft - ft 1 *) 2 + 7 - n) 2 ] 

i=i (Tu ' / i = 1 

for the model y = /3q + /3\x. This expression assumes that the variance ratio is 
a known value. Any errors in estimating the variance ratio are neglected. Using 
standard calculus methods, the estimates of the model parameters are 


Po = Y - PiW, 


Pi = 


S Y y — 7 Sww + \] (S Y y — 7 Sww ) 2 + 47*5 % Y 


25, 


WY 


, and 


Xi = Wi + 


§1 

Pi + If 


Yi-Po-PiWi 
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where 


Y 


1 

n 




W 


1 

n 




Syy = J2{Y<~ Y) 2 , 
i = 1 


Sww = ^2(Wi-W) 2 , 

i = 1 

and 

n 

S W Y = ]T (w t - W) (Yi - Y) . 

i — 1 

When 7 = 1 , the estimates of the model coefficients are equivalent to the OrthLS 
estimates. The only difference is that when W is unknown, x is substituted into these 
expressions. It is important to note that for normally distributed errors, the Deming 
estimates are also maximum likelihood estimates (Casella and Berger, 2002). 

Deming regression is considered a special case of total least squares (TLS) (Golub 
and Loan, 1980), which is used for multivariate ME models. Unlike Deming regres- 
sion, TLS does not generally have a closed-form solution. The method is a numerical 
solution to the linear algebraic equation 


(x + u) f3 = y + e 


where the objective is to minimize u and e given x and y. Singular-value decompo- 
sition is one approach to solving for the model coefficients in the equation. However, 
implementing TLS in practical application, such as calibration, is difficult. 

Carroll and Spiegelman (1986) discuss the effect of small MEs on the estimates 
the regression coefficients and confidence intervals in an instrument characterization 
experiment. In their discussion, Carroll and Spiegelman only consider the simple 
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linear model. For most calibration applications at NASA LaRC, the MEs are at least 
four times smaller than the accuracy of the measurement system being characterized. 
Carroll and Spiegelman acknowledged two criterion used to determine whether ME, 
or uncertainties in the standard, can be neglected in the analysis. These criterion 
are: 


• the ratio of ME to the true variability in the x’s (Draper and Smith, 1998), 
and 

• the ratio of the measured response to the slope, or primary sensitivity (Scheffe, 
1973; Mandel, 1984). 

They stated that the first criterion is most appropriate for estimating the regression 
coefficients while the second criterion affects the width of the confidence intervals, 
but they argue both criterion should be used in combination when analyzing the 
data. 

The ratio, u, based on Draper and Smith’s (1998) criterion is 

ME Variance o\ 

Variance of x S xx 

In the limit as the number of observations, n, approaches infinity, the estimated 
intercept is approximately 

A> ~ A + 

and the estimated slope is approximately 

h ~ (At) 


Therefore, for small u, the estimated coefficients are only affected slightly by ME. 
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Centering the x’s also negates any effect of ME on the estimate of the intercept. 
The width of the confidence interval without ME is approximately 2t a /2, n -p^e, 
where t a / 2 , n - P is the (1 — a/2) percentage point of the ^-distribution with n — p 
degrees of freedom. When MEs are ignored, the confidence interval is approximately 

1 I r\ 

2t a /2,n- P [pi + Pi&u) . As a result, the ratio of confidence interval width when 
ignoring ME to conhdence interval width without ME is 

2 t a /2,n-p (of + P\ fap 1/2 = f 1+ * 

2^a/2,n— p&e \ ®1 / 


p 2 a 2 \ 1/2 

It is evident that the conhdence interval width is increased by I 1 H — | when 


ignoring the ME. This result is also true for prediction intervals. 

Carroll and Ruppert (1996) study the effect of modeling error in OrthLS. In any 
regression, it is desired that the estimated mathematical model adequately char- 
acterize the relationship in the experimental data. Standard regression techniques 
utilize a lack-of-ht (LOF) test to assess model adequacy, but this test is a function 
of the pure experimental error and the available number of degrees of freedom. In 
experiments with MEs, replicates are not genuine; instead, data points that are close 
enough may be considered pseudo- replicates. Therefore, without a LOF-equi valent 
test for ME models, it is difficult to assess model adequacy. Carroll and Ruppert go 
on further to suggest that incorrect knowledge of the variance ratio also attributes to 
LOF in the model. They argue that in practice 7 is underestimated, which leads to 
an overcorrection in the estimated slope in OrthLS. For appropriate use of OrthLS, 
the variance ratio should be redefined to be 


7 = 


^ 2 + ffi 2 
cC 



21 


where a ^ is the variance of the modeling error and is assumed to be constant. In simu- 
lation, the modeling error may be neglected since the appropriate model form is often 
known. For practical applications, Carroll and Ruppert recommended that OrthLS 
be used with caution, especially when the appropriate model form is unknown. 
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CHAPTER 3 

SIMPLE LINEAR MODELS 

3.1 APPLICATION OF THE SIMPLE LINEAR MODEL 

In the literature review, a body of research associated with measurement errors 
(ME) for the simple linear model was discussed. For some applications, the simple 
linear model works well. In general, the simplest model that adequately explains the 
experimental data is desired. It is important to recognize that any estimated model is 
only an approximation of the true, underlying function. There are many applications 
that require multi-dimensional, higher-order response surface models, and the theory 
for the simple linear problem provides the necessary framework on which to expand 
in order to model more complex situations. Therefore, the discussion begins with 
the simple linear model in order to formulate a general methodology. 

In this research, some of the most common techniques used for MEs in the simple 
linear model were explored. While additional statistical theory is available, the prac- 
tical implementation of this theory is limited. As a result, the goal of this work was 
to develop an engineering solution using some principles from least squares theory. 
The results are applicable to practical situations including the calibration of simple 
instruments such as single-axis load cells and pressure transducers. 

Consider an instrument where the output, y, is assumed to be linearly related to a 
single applied factor, x. A calibration experiment is designed to measure the response 
at known values of x. The resulting set of n design points are used to estimate the 
linear relationship between x and the observed response Y. The estimated model is 
then employed on future observations of x to predict Y, as long as the value of x falls 
within the region spanned by the calibration experiment. However, in the calibration, 
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it is assumed that the applied values of x are known, fixed quantities. In practice, 
the actual value of x applied is seldom known without error. Since OLS estimation 
requires that the regression parameters be known without error, it is necessary to 
utilize practical assumptions with regard to their actual accuracies. The two most 
common assumptions are: 

• The error in x is small relative to the error in the response. In calibration 
applications, the magnitude of the error in x is typically stipulated. 

• The error in x is captured in the error in the response. This results in the 
model containing a single error term; namely e* = e + u. 

For ME models, Fuller (1987) demonstrated the consequences of this naive analysis, 
including induced bias in the estimate of the slope coefficient. 

As a motivating, practical example for the simple linear case, consider the NASA 
Orion Multi-Purpose Crew Vehicle (MPCV), which is currently being designed to 
transport astronauts beyond low-Earth orbit (LEO). The first flight test of the Orion 
MPCV will be the Exploration Flight Test 1 (EFT-1) and will be a four-hour orbital 
test scheduled for late 2014. The flight test will test various spacecraft subsystems 
and includes several on-board measurement systems to collect flight data. One of 
these measurement systems is a Flush Air Data System (FADS). The FADS for 
EFT-1 contains nine pressure ports located in a cruciform pattern across the fore- 
body heatshield of the aeroshcll and the data collected will be inputs to post-flight 
trajectory reconstruction algorithms that estimate various flight parameters, such 
as angle-of-attack and Mach number. Each pressure port is connected via tubing 
to a digital pressure transducer that is mounted on the backside of the heatshield. 
The uncertainty in estimating the flight parameters is a function of the calibration 
of each pressure transducer (Commo and Parker, 2012). The transducer output, y , 
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is assumed to be linearly related to the sensed pressure, x. To verify this assump- 
tion, a simple calibration experiment was designed to estimate the mathematical 
relationship between the output and sensed pressure. The pressure levels during the 
calibration were applied via a reference pressure standard. Reference standards at 
NASA are calibrated against a NIST-traceable standard. However, reference stan- 
dards contain a small amount of uncertainty although they are often used in a similar 
manner as a NIST standard. As a result, the applied pressures are only known to 
within a certain accuracy. The actual applied pressure in this case is 

W = x + u 

where u is the uncertainty in the reference standard as estimated from the NIST- 
traceable calibration. 

Assuming the simple linear model is appropriate for this example, then 


Vi = Po + Pi Xi + 6i (3) 

where p 0 is the zero-offset, Pi is the sensitivity, and e is the normally-distributed 
random error of y with a mean of zero and a constant variance of af . The stan- 
dard deviation, a e , is typically considered to be the accuracy of the instrument in 
the classical model, but the accuracy is redefined later for the ME model. As ex- 
pressed, Equation (3) is consistent with the ads being known without error. Under 
this assumption, OLS is appropriate for estimating the regression coefficients. In the 
presence of MEs, Equation (3) becomes 


Vi = Po + Pi {Xi + Ui ) + 6i 


( 4 ) 


where Ui is the error in ay and is distributed with a mean of zero and a constant 



25 


variance of a 2 u . If u was known for every instance of x , then estimation could proceed 
with OLS estimation using W = x + u as the regressor instead of x. In physical 
experiments, the value of u is not known and other methods are required to estimate 
the regression coefficients. 

From statistical regression theory, OLS minimizes the squared difference between 
the observed value of the response, Y, and the predicted observed value of the re- 
sponse from the estimated model, Y . Mathematically, this is expressed as 


An $1 


lb 2 

min ^2 (A - • 

i — 1 


If the observed value contains no error, then the minimization is 


An A 


n 

min E (Vi-Vi) 2 ■ 

i = 1 


It is easily seen that any MEs are not considered in this objective function, therefore, 
the objective function is one-dimensional. For the simple linear model with MEs, a 
two-dimensional objective function is appropriate. The form of this objective function 
is determined by the variance ratio, 7. The variance ratio is 


<7; 


7 = 


at 


This ratio is particularly important in determining how well the ads should be known 
in a calibration experiment. For OLS, 7 is assumed to be zero. 

One particularly well-known and useful estimation technique in the presence of 
ME is orthogonal least squares (OrthLS) (Carroll et ah, 1985). Most of the literature 
available on OrthLS is related to the simple linear problem. While OrthLS provides 
a nice alternative to OLS for ME models, it relics on the assumption of equal error 
variances of both variables. In other words, the variance ratio is assumed to be one. 
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In terms of the calibration application, the uncertainty in the calibration reference is 
on the order of the uncertainty in the response, which is insufficient. It is preferred 
that the uncertainty in the reference standard be at least four (4) times smaller than 
the uncertainty in the response of the measurement system. This corresponds to a 
variance ratio of 



The variance ratio for typical calibration applications at NASA are 0.0625. Because 
of this, neither OLS or OrthLS is most appropriate for estimation of the mathematical 
model. However, there is no known least squares-based technique that can be applied 
for variance ratios between 0 and 1. 

The uncertainties in the responses of the measurement systems at NASA can vary 
drastically. Highly- accurate instruments, such as the pressure transducers for EFT-1, 
range from 0.01 to 0.1 percent of the full-scale (F.S.) range. For a 0.01 percent F.S. 
uncertainty, the calibration reference should have a maximum uncertainty of 0.0025 
percent of F.S. Other measurement systems, like force-balances, are transducer-class 
instruments with slightly larger uncertainties, on the order of 0.10 to 10 percent of 
F.S. The worse-case is seen in atmospheric science applications where uncertainties 
can approach 50 percent of F.S. For this research, the focus is limited to a maximum 
uncertainty of 20 percent of F.S. 

For the simple linear model, Figures 1-4 show the spatial dispersion of errors in 
x and y for several variance ratios. It is obvious that as errors in x become larger, 
the assumptions required for OLS estimation no longer hold and the estimate of the 
slope is affected. Therefore, other estimators, such as OrthLS for 7 = 1, are more 
appropriate and more robust with regards to the MEs. In addition, the implications 
of the two ME assumptions discussed earlier are revealed. For small MEs relative to 
the errors in the response (i.e. 7 = 0.01), Figure 2 shows that the estimate of the 
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slope can still be affected. If the ME is combined with the response error to form 
a single error term, then the dimension of error is reduced from two to one, as in 
Figure 1. 
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FIG. 1: Spatial Dispersion of Errors for the Simple Linear Model for 7 = 0 



FIG. 2: Spatial Dispersion of Errors for the Simple Linear Model for 7 = 0.01 
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FIG. 3: Spatial Dispersion of Errors for the Simple Linear Model for 7 = 0.25 



3.2 ESTIMATION METHODS FOR THE SIMPLE LINEAR MODEL 


For the simple linear measurement system characterization problem, the following 
three methods of estimating the model coefficients are discussed: OLS, OrtliLS, and 
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propose modified least squares (MLS) as a more general approach to address the 
limitation previously mentioned. 

3.2.1 ORDINARY LEAST SQUARES 


For a calibration experiment with n design points, the sample means are 


and 


x = ~y ^ Xi 
1=1 


r = ±I> 

n 


i = 1 


where the responses are the error-prone values. Using the sample means, the sum of 
squares are 


S X x ^ ^ {Xi 


X) 


1=1 


and 


s„ = E( y <- ? ) 2 - 
1=1 

Similarly, the sum of cross-products is 

n 

s xy =yy j {.x i - x) (Yi - y) . 

2=1 

As stated previously, the minimization criteria for OLS is 

n 

'EM- Po-Pi 




fjL / . N 2l 


An Al 

= min 

K 

1 

= min 



Li=i ' ' J 



Xj 


i= 1 


(5) 


The OLS estimators are found by differentiating the objective function in Equation 
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(5) with respect to A) and Al and setting each resulting equation to zero. Therefore, 

- o i n 

w =- 2 E( r <-A-M=° 

i= i v 7 

and 

"(9 1 n 

These are known as the least squares normal equations. Solving for (3 0 and Al, the 
OLS estimators for the simple linear model are 

A) = Y - fax (6) 

and 

n 

E (x i -x)(Y i -Y) 

h = ^ ( 7 ) 

E (Xi - xf 

X=1 

When the xs are centered, x = 0 and the estimator of the zero- intercept, /3q, is Y . 
As mentioned earlier, the Gauss-Markov Theorem states that under the objective 
function given in Equation (5) and for uncorrelated, homoscedastic errors, the esti- 
mates /3o and Al are best linear unbiased estimates (BLUE). Simply put, under the 
criteria of minimizing the distance between Y and Y, the OLS estimates will always 
yield the smallest sum of squared errors. This concept leads to a new definition of 
residual error for ME models as a function of the variance ratio. 

3.2.2 ORTHOGONAL LEAST SQUARES 

When the xs are known without error, the OLS criteria for minimizing the errors 
in the ys is a logical choice. However, in the ME models, it is more appropriate to 
use an estimator that considers the error in the xs in addition to the ys. The most 
commonly-employed method that accounts for errors in both variables for the simple 
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linear model is OrthLS. In OrthLS, the estimators are found by minimizing the sum 
of the squared orthogonal distances from all the points to the estimated line. See 
Figure 5 for the geometrical interpretation of the two estimators. Therefore, OrthLS 
equally weighs the errors in both x and y when deriving the estimators. In terms of 
the variance ratio, 7 = 1 . 



FIG. 5: Distance Minimized by Orthogonal Least Squares 


Now consider the data point (aq,li) in Figure 5. Since u is unknown and an 
additive ME model is assumed, the best estimate of the unknown, actual level of 
factor W is x. Therefore, the derivation of the OrthLS estimators proceeds based on 
x. From algebra, it is shown that the shortest distance from a point to a line is the 
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orthogonal distance. If the line is given by the simple linear model 

Y — A + PiX, 

then the corresponding closest point on the line is (Stewart, 2003) 

- A A + %i — AA , Q s 

x ' = TWi ’ (8) 

and 

Y\ = A + ~7~~R2 (A^ i + x i ~~ AA) ■ (9) 

i + Pi 

The squared orthogonal distance between the point and the line is defined as 
(x — x) 2 + — Y j . Over n design points in an experiment, the objective function 

is 

n r 2" 

A, A = min X] + (*< “ Yij ■ 

i=l ^ 

Substituting Equations (8) and (9) into the objective function 

, . T 1 n 

A, A = min ! , 02 X [A _ (A) + A®*)] 2 • (io) 

L J [ i+Pl i=l 

Differentiating Equation (10) with respect to A) and setting the resulting equation 
to zero yields A — Y — (3\X, which is identical to the OLS estimate of the zero- 
intercept. Regardless of whether the errors in x are included or excluded, the estimate 
of the zero-intercept is unaffected in the simple linear model. Substituting this result 
back into the right-hand side of Equation (10) and setting equal to zero to find the 
minimum gives 

-t Yi 2 1 

^ X \( Yi - “ A ( x i - ®)1 = j ^ \ s yy ~ 2 A s xy + AiAJ = 0 

1 + A ~i 1 J i + A L J 
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using the prior definitions of the sample sum of squares and sum of cross-products. 
The minimum and thus the OrthLS estimate of the sensitivity is 


A 



s t 


yy ) 


+ 


( 11 ) 


which is similar to Deming’s estimate of the slope. As in the OLS derivation, we 
define the residual error as the distance between a point and the estimated line. 
Hence, the residual error for OrthLS is (x — x) 2 + — Y ^ . 

3.2.3 MODIFIED LEAST SQUARES 

The estimators for two cases, 7 = 0 and 7 = 1 , have been derived. As men- 
tioned previously, for most applications, the variance ratio is between 0 and 1. How- 
ever, a least squares estimator that can be used generally for these cases would be 
useful. The geometry of the problem is exploited to help define the new estima- 


tor. 
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FIG. 6: Distance Minimized by Modified Least Squares 


Consider the point (aq,Yi) in Figure 6. The vertical line segment between the 
point and the line is a = Y\ — Y*. The angle, </>, is made by the line segment of length 
a and the estimated line as 


0 = 


n _i 

2 “ tan ~ 


dy 

dx 


dy 

where — is the slope of the line or f3\. Next, the line segment d is defined as the 
dx 

distance minimized in MLS. This distance is a function of both the slope of the 
estimated line and the variance ratio. From d and the estimated line, a second angle, 
a, is defined as 

<* = ! + (! -7)tan-T| 


where 7 is the variance ratio. The angle a is a first-order approximation based on 
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known angles when 7 = 0 and 7 = 1. From trigonometry, the Law of Sines states 

sin a sin <f> 
a d 


Solving for d yields 

sin 0 

d = a. 

sm a 


( 12 ) 


Using Equation (12), the newly defined objective function for MLS minimizes the 
sum of the squared lengths of the dd s or 


An Pi 


min di = min 


i— 1 


Z_. 

i= 1 


sm 0 

-( 

sina 


n 2 


(13) 


Note that as a — > |, Equation (13) is equivalent to the OrthLS estimator. As in the 
case of OrthLS, the only additional information required to use Equation (13) is an 
assumption about the variance ratio. For the calibration applications, the variance 
ratio varies between 0.0625 and 0.0001. 


3.3 SIMULATION STUDY FOR THE SIMPLE LINEAR MODEL 

The statistical properties of the OLS estimators are well-known. From the Gauss- 
Markov theorem, the expected values of Pq and Pi are A) an d Pi, respectively due to 
the unbiased nature of the OLS estimator (Myers, 1990). Furthermore, the variance 
of the estimators are 

Var (*) = aj (i + £) 

and 

2 

Var (A) = 

^ ' ^xx 


and these are minimum variance estimators. To understand the statistical properties 
of the OrthLS and MLS estimators, a simulation study was conducted to make 



36 


inferences on the effects of MEs on the prediction capabilities of each method. The 
prediction capabilities are based on the residual error, which has been redefined to be 
the variance ratio-weighted distance between any point and the estimated line. This 
is represented by the distance, d, in Figure 6. The simulation study was designed 
to be representative of an actual calibration used on simple instruments, such as 
the pressure transducers for EFT-1 discussed earlier. A calibration experiment is 
designed with n = 6 design points in order to fit the simple linear model 


y = (3 0 + (3 1 x + e. 


In the presence of ME, this model becomes 

y = p 0 + (3i (x + u) + e. (14) 

The following assumptions were made with regard to Equation (14): 

• /3 0 — 0 and (3\ — 1. Through appropriate centering and scaling of the data, 
this is always a plausible option in real experiments. 

• The errors u and e are independent and identically distributed as normal with 
means of zero and constant variances of of and of, respectively. 

• of and of are proportionally related through the variance ratio, 7. 

• Wi and Y{ are the error-prone values of 07 and ?/*. Therefore, Wi and Y t are 
jointly distributed as a bivariate normal distribution 


[W i ,Y i \~N[(x i ,(3 Q + f3 1 x i ),Y] 


where E = diag (yof , of). Both Wi and Y t are independent random variables. 
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From these n design points and corresponding responses, the zero- intercept and sen- 
sitivity coefficients were estimated, and the mean squared residual errors (MSE) were 
calculated for each of the estimation methods. The MSE is defined as 

MSE = (Bias) 2 + Precision 

Next, five (5) confirmation points located at x — —1, —0.5, 0, 0.5, and 1 were used 
to generate confirmation data to further examine the prediction capabilities of each 
estimated model. For a given calibration design, each method estimates a set of 
coefficients M times. Each set of estimated coefficients are used L times on the 
confirmation points and the MSE is calculated. 

The simulation study considered three (3) different calibration designs as shown 
in Table 2. Each design is sufficient for estimating the regression coefficients in the 
model, and each design has certain advantages and disadvantages over the other de- 
signs. Design #1 is a balanced design with replicates throughout the design space. 


TABLE 2: Designs Considered for the Simple Linear Simulation Study 


Design 

Design Points 

Comments 

Design #1 
Design #2 
Design #3 

-1,-1, 0, 0, 1,1 
-1, 0, 0, 0, 0, 1 
-1,-1, -1, 1, 1, 1 

Replicates throughout Design Space 
Suspicion of Significant Lack-of-fit 
Strictly First-order Design 


This allowed testing of the constant variance assumption, and for calibration appli- 
cations, this replication strategy provided insight on the repeatability of the mea- 
surement system over the calibration range. Design ^2 contains the same number 
of degrees of freedom for replication as Design #1, but it does not allow testing of 
constant variance throughout the design space. This design is suited for instances 
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when lack-of-fit is suspected to be significant. Both designs allow for estimation of a 
second-order model term. Design ^3 is a first-order-only calibration design. Repli- 
cation is dispersed equally at the low and high ends of the range so testing of the 
constant variance assumption is permissible. However, inferences on the estimation 
methods are not design-dependent. Table 3 shows the values of 7 and a e that are used 
as inputs to the simulation. For variance ratios, the maximum variance ratio stud- 
ied is 7 = 1, representing errors of equal orders of magnitude. Commonly-assumed 


TABLE 3 : Variance Ratio and Response Uncertainties Considered for the Simple 
Linear Simulation Study 


Run No. 

Variance Ratio, 7 

Response Uncertainty, cr € , % F.S. 

1 

1 (1:1) 

0.2 ( 20 % F.S.) 

2 

0.25 ( 2 : 1 ) 

0.2 

3 

0.0625 ( 4 : 1 ) 

0.2 

4 

0.01 (10:1) 

0.2 

5 

0.0001 (100:1) 

0.2 

6 

1 

0.1 ( 10 % F.S.) 

7 

0.25 

0.1 

8 

0.0625 

0.1 

9 

0.01 

0.1 

10 

0.0001 

0.1 

11 

1 

0.05 ( 5 % F.S.) 

12 

0.25 

0.05 

13 

0.0625 

0.05 

14 

0.01 

0.05 

15 

0.0001 

0.05 

16 

1 

0.01 ( 1 % F.S.) 

17 

0.25 

0.01 

18 

0.0625 

0.01 

19 

0.01 

0.01 

20 

0.0001 

0.01 

21 

1 

0.001 ( 0 . 1 % F.S.) 

22 

0.25 

0.001 

23 

0.0625 

0.001 

24 

0.01 

0.001 

25 

0.0001 

0.001 
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ratios for most calibration applications at NASA are also included in the simulation, 
such at 7 = 0.0001, 0.01, 0.0625, and 0.25. The response uncertainties studied in this 
simulation are representative of a large population of measurement systems at NASA 
so the F.S. accuracies range from 0.1 to 20 percent. The simulation study employed 
a factorial combination of variance ratios and response uncertainties. An outline of 
the simulation study is provided below: 

1. For a given calibration design in Table 2, generate the yj s from the simple 
linear model /3 0 + PiXi, where /3 0 = 0 and Pi = 1. 

2. Using ( Xi,yi ) for i = 1,2, ..., 6, generate (Wj,Y)) from a bivariate normal 
distribution. The distribution is [1'U, Y t ] ~ N [(ay, Po + PiXp , E], where E = 
diag (yof, of) and 7 and oy are taken from Table 3. 

3. From (ay, Y)) for i = 1, 2 , . . . , 6, estimate ^/3o, Pi j for OLS, MLS, and OrthLS. 

Using compute the predicted value y r for each method and calculate 

the MSE based on 7 and oy. 

4. For the confirmation points x c = [—1.0, —0.5, 0.0, 0.5, 1.0], generate the y Ci ' s 
from the simple linear model Po + P\X Ci , where Po = 0 and Pi = 1. 

5. Using ( Xi,yp for i = 1,2,..., 5, generate ( W Ci ,Y Ci ) from a bivariate normal 
distribution. The distribution is [IF^, Yf ] N [(xi, Po + PixP , E], where E = 
diag (7 o’g , al) and 7 and oy are the same values as in Step 2. 

6. Using compute the predicted value Y Ci for each method and calculate 

the MSE based on the values 7 and oy in Step 2. 

7. Repeat Steps 5 and 6 for L = 1000 times. 

8. Repeat Steps 2 through 7 for M = 100 times. 


9. Repeat Step 8 for K = 25 times. 
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3.3.1 SIMULATION RESULTS 

Tables 4 shows the mean estimate of the sensitivity coefficient for the three cali- 
bration designs and the three estimation methods. The mean estimates of the zero- 
intercept are not provided since all three methods yield the same value for a given 
design, variance ratio, and response uncertainty. From Table 4, the following infer- 
ences are made based on the structure of the simulation: 

• For a variance ratio of one (7 = 1), the estimates of j3\ obtained from MLS and 
OrthLS are equivalent as expected. In the limit as 7 — > 0, the estimate of f3\ 
from MLS converges towards the OLS estimate. 

• For an uncertainty of 1 percent or less < 0.01), the three methods yield 
nnmerically-similar estimates of the sensitivity coefficient to four significant 
figures. 

• Larger differences between the estimates are observed when the uncertainty 
is worse than 5 percent ( a € > 0.05). The magnitude of the difference varies 
between designs and levels of response uncertainty. 

• The difference between the OLS estimate and OrthLS is approximately constant 
across the levels of variance ratio for a given design and uncertainty. For 
example, in the first design with an uncertainty of 20 percent, the differences 
are 

7 = 1 : A = 1.0229 - 1.0031 = 0.0198, 

7 = 0.25 : A = 1.0177 - 0.9975 = 0.0202, 

7 = 0.0625 : A = 1.0172 - 0.9967 = 0.0205, 

7 = 0.01 : A = 1.0173 - 0.9974 = 0.0199, and 


7 = 0.0001 : A = 1.0231 - 1.0033 = 0.0198. 
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In this example, the differences correspond to a range of approximately 2 per- 
cent. 

Table 5 presents the variance of the mean estimate of the sensitivity coefficient for the 
three designs and the three methods. For a given design and level of uncertainty, the 
variance of the estimates were consistent between the methods and across the levels 
of variance ratio. Small differences were observed in the magnitude of the variance 
between the designs for any level of uncertainty, but this could be attributed to the 
design properties rather than the estimation method. However, from a practical 
application perspective, using MLS to estimate the sensitivity coefficient resulted in 
the same magnitude of variability as OLS. 



TABLE 4: Mean Estimate of the Sensitivity Coefficient for the Simple Linear Simulation 
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TABLE 5: Variance of the Mean Estimate of the Sensitivity Coefficient for the Simple Linear Simulation 
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Table 6 shows the average MSE across all the design and confirmation points for 
the three designs and the three methods. Recall that the MSE includes both bias in 
the prediction and precision of the prediction. The root mean squared error (RMSE) 
is the square root of MSE and is typically used to quote the estimated accuracy of a 
measurement system after calibration. For example, the 1-er estimated measurement 
system accuracy for Design #1, 7 = 1, and a e = 0.2 based on OLS is 


/ Mean Squared Error V0.0132 
Full-Scale Range 1 


or 11.5 percent of F.S. While Table 6 shows the absolute MSE values, Table 7 reveals 
the difference between OLS and either MLS or OrthLS. From this table, it is evident 
that MLS has a smaller MSE than OLS in all the cases studied. Additionally, OrthLS 
shows a benefit over OLS in most cases when the variance ratio is greater than or 
equal to 0.25. These differences are expressed in terms of percent improvement in 
Table 8. Both MLS and OrthLS provide up to a 3 percent reduction in MSE over 
OLSE, such as in the case of 7 = 1 and cr e = 0.2. As the variance ratio departed 
from unity, the relative improvement of MLS decreased while OrthLS showed neg- 
ative impact at small values of 7. Through the use of one-way analysis of variance 
(ANOVA), the differences observed between the methods were tested for statistical 
significance. ANOVA compares the means of two or more groups in the presence of 
the variability within each of the groups (Casella and Berger, 2002). Using a criteria 
of 0.05, the hypothesis that the MSE between methods are equal is rejected if the 
p-value from the ANOVA is less than 0.05. ANOVA does not directly indicate which 
groups are different or by how much the groups differ. The ANOVA is shown in Ta- 
bles 9-11. Separately, the tables reveal that the difference in MSE is not statistically 
detectable since the p-value is near one. When comparing the three hypothesis tests 
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simultaneously, the rejection criteria becomes 

atotai = 1 - (1 - « indiv ) n = 1 - (1 - 0.05) 3 = 0.14 

where n is the number of simultaneous tests. This value, also known as the family- 
wise error rate, is the probability of making a Type 1 error across all the hypothesis 
tests. However, the conclusion remains that there is no detectable differences in the 
MSE between the three methods. 

Table 12 emphasizes a practical interpretation of the differences in the MSE. The 
percent improvement in MSE relative to the uncertainty in the response is calculated 
by 

^Improvement in Mean Squared Error 

If the improvement in MSE is less than zero, then the percent improvement is zero. 
These results show that the impact of MLS is greater when comparing the improve- 
ment to the uncertainty in the response for the three calibration designs studied. For 
example, a 0.01 percent improvement in MSE is achieved on a measurement system 
with a 0.1 percent response uncertainty using the first calibration design and a vari- 
ance ratio of 0.0625. This difference is considered practically significant. In summary, 
while the statistical hypothesis tests revealed no difference in the MSE between the 
methods, the practical benefit of MLS is evident in Table 12. For variance ratios 
greater than 0.0625, it is recommended that MLS be employed instead of OLS. 



TABLE 6: Mean Squared Error for the Simple Linear Simulation 
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TABLE 7: Improvement in Mean Squared Error over OLS for the Simple Linear Simulation 
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Note: Asterisk represents an improvement of less than 0.00005 




TABLE 8: Percent Improvement in Mean Squared Error over OLS for the Simple Linear Simulation 
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Note: Asterisk represents an improvement of less than 0.005% 
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TABLE 9: Analysis of Variance of the Mean Squared Error for the Simple Linear 
Simulation - Design #1 


Source 

DF 

SS 

MS 

F 

P 

Method (OLS, MLS, OrthLS) 
Error 

2 

447 

0.0000002 

0.0358240 

0.0000001 

0.0000801 

0.00 

0.999 

Total 

449 

0.0358242 





TABLE 10: Analysis of Variance of the Mean Squared Error for the Simple Linear 
Simulation - Design ^2 


Source 

DF 

SS 

MS 

F 

P 

Method (OLS, MLS, OrthLS) 
Error 

2 

447 

0.0000008 

0.0397981 

0.0000004 

0.0000890 

0.00 

0.996 

Total 

449 

0.0397989 





TABLE 11: Analysis of Variance of the Mean Squared Error for the Simple Linear 
Simulation - Design ^3 


Source 

DF 

SS 

MS 

F 

P 

Method (OLS, MLS, OrthLS) 
Error 

2 

447 

0.0000001 

0.0342407 

0.0000000 

0.0000766 

0.00 

0.999 

Total 

449 

0.0342408 









TABLE 12 : Percent Improvement in Mean Squared Error Relative to Response Uncertainty for the Simple Linear Simulation 
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Note: Asterisk represents an improvement of less than 0.005% 
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FIG. 7: Percent Improvement in Mean Squared Error of MLS for the Simple Linear 
Simulation - Design #1 
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FIG. 9: Percent Improvement in Mean Squared Error of MLS for the Simple Linear 
Simulation - Design #3 
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CHAPTER 4 

MULTI-DIMENSIONAL, HIGHER-ORDER MODELS 

4.1 APPLICATION OF THE SIMPLE POLYNOMIAL MODEL 

In the previous chapter, the simple linear model was discussed in detail, and mod- 
ified least squares (MLS) was proposed as an alternative method to ordinary least 
squares (OLS) for estimating the regression coefficients in the presence of measure- 
ment errors. The statistical properties of the estimator and the prediction capabilities 
of MLS were studied through simulation along with OLS and orthogonal least squares 
(OrthLS). Results showed that the numerical differences in estimating the sensitiv- 
ity coefficient were small between the three methods. Furthermore, differences in 
the mean squared error (MSE) were also small between the methods, and through 
the use of analysis of variance (ANOVA), the differences were deemed statistically 
undetectable. However, the practical implications of the different methods were pro- 
nounced. For the three designs considered, MLS showed a substantial improvement 
in MSE over OLS relative to the uncertainty in the response. As an example, MLS 
had a 0.06 to 0.51 percent reduction for an response uncertainty of 0.01, depending 
on the level of the variance ratio. Based on these results, it was recommended to 
employ MLS for simple linear models when the variance ratio was larger than 0.0625. 

While the simplest form of the estimated model is desired, many applications are 
more complex and include additional factors or higher-order model terms. For now, 
the focus is limited to single-factor, polynomial models. The multiple-factor case 
is discussed shortly. The single-factor, polynomial model is identified as the simple 
polynomial model within this research. For OLS estimation of a simple polynomial 
model, the theory is readily available and is briefly reviewed in this chapter (Myers, 
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1990). Methods comparable to OrthLS for simple polynomial models are under- 
developed within the literature. Other estimation techniques, such as those based 
on the Method of Moments, for the simple polynomial model are limited in practical 
application and available in popular ME references (Fuller, 1987; Buonaccorsi, 2010). 
The methods of OrthLS and MLS discussed in the previous chapter provide the 
necessary framework to expand upon for simple polynomial models. Developing 
these methods for other model forms is important to understanding the implications 
of using naive estimation methods, such as OLS. 

In this research, the focus is on second-order polynomial models, which are used 
in several calibration applications. Similar to the simple linear problem, the goal is 
to develop an engineering solution that can be easily applied to practical situations 
where the estimated model is a polynomial. Consider again the calibration of the 
pressure transducers for EFT-1. Due to certain properties of the sensing element 
inside the transducer, the applied pressure is suspected to have a nonlinear effect on 
the response. The magnitude of this nonlinear effect is small relative to the primary 
sensitivity of the transducer but nonetheless is significant in affecting the response. 
Therefore, the simple linear model may not be sufficient in capturing the functional 
behavior between x and y. A set of n design points are used to estimate the regression 
coefficients in the model, and the assumed appropriate form of the model is now 

Vi — Po + fiiXi + /3 uX} + (15) 

where /3 0 is the zero-offset, /A is the sensitivity, (i\ \ is the nonlinear effect, and e is the 
normally-distributed random error of y with a mean of zero and a constant variance 
of of. The applied pressures are only known to a certain accuracy and as a result, 
the actual applied pressure is 


W = x + u 
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where u is the uncertainty in the reference standard as estimated from the NIST- 
traceable calibration. With MEs, Equation (15) becomes 


Vi = Ao + Al + Ui) + All (Xi + Uif + €i (16) 

where u t is distributed with a mean of zero and constant variance of Since the 
value of u at a value of x is unknown, the variance ratio is used to develop both MLS 
and OrthLS as alternative estimation methods to OLS in the presence of MEs. 

When comparing the effect of MEs on the simple linear and quadratic models, 
there are a few important points to remember. While the assumed form of the 
ME model is additive, the expected value and variance of the response are affected 
by the order of the model, as discussed earlier in Chapter 1. For convenience, the 
results are summarized for the simple linear and quadratic models below in Table 
13. In the simple linear model, the expected value for classical and measurement 
errors are similar while the variance is inflated due to the ME. Conversely, both the 
expected value and variance of y are affected by the ME in the quadratic model. The 
expected value of y is biased by AiM and the variance of y is further inflated by 

oi( 1 + 4/3M). 


TABLE 13: Comparison of Bias and Variance of Simple Linear and Quadratic Models 




Classical Error 

Measurement Error 

Simple Linear Model 

Ely) 

Ao + Pixi 

Ao + Ai®i 

Var (y) 


of + AM 

Simple Quadratic Model 

Ely) 

A) + Ai®i + An®? 

Ao + Ai®i + An ( x i + a u) 

Var (y) 

of 

+ Of (1 + A? + 4AM?) 
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The objective function within each method is a geometric interpretation. For 
instance, OLS minimizes the sum of the squared differences between the observed 
response, Y t , and the predicted response, Y . over n design points. Geometrically, this 
is equivalent to the E-distance between a point and the estimated line. Alternatively, 
OrthLS minimized the sum of the squared orthogonal distances between n design 
points and the estimated line. MLS minimizes the sum of the squared variance- 
weighted distances and is bound at the extremes by OLS and OrthLS. One advantage 
of OLS in practical applications is that the minimization yields a closed-form set 
of equations that are readily solvable. The other two methods lend themselves to 
numerical solutions. Because of this, OLS is simpler to use in practice and most 
basic software packages contain OLS regression. 

For the simple quadratic model given by Equation (15), Figures 10-13 demon- 
strate the effect of variance ratio on the distribution of errors in x and Y . In Figure 
10 where 7 = 0, it is obvious that the errors are again one- dimensional and because 
of this, OLS is a sufficient estimator of the model. Additionally, while the curvature 
is small, it is still observed in Figure 10. As 7 — > 1 , it is seen that both the slope and 
curvature are affected by the errors in x and OLS is no longer a sufficient estimator 
since the errors are two-dimensional. In the simple linear case, only the slope was 
affected by any MEs. Furthermore, the errors in x can also disguise any departures 
from linearity, as seen in Figure 13. 
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FIG. 10: Spatial Dispersion of Errors for the Simple Quadratic Relationship for 7 = 0 



x 


FIG. 11: Spatial Dispersion of Errors for the Simple Quadratic Relationship for 
7 = 0.01 
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FIG. 12: Spatial Dispersion of Errors for the Simple Quadratic Relationship for 
7 = 0.25 



FIG. 13: Spatial Dispersion of Errors for the Simple Quadratic Relationship for 7 = 1 


4.2 ESTIMATION METHODS FOR THE SIMPLE QUADRATIC 
MODEL 

For the simple quadratic measurement system characterization problem, the 
methods of OLS, MLS, and OrthLS are again discussed. As a note, the OLS estima- 
tor derived shortly is applicable to both polynomial and multi- dimensional response 
surface models. 


4.2.1 ORDINARY LEAST SQUARES 
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Consider the calibration of a measurement system with k factors of interest and n 
design points are executed during the calibration. These n design points are sufficient 
in estimating a model of the general form 

Yi = (3 [ o + + . . . + / 3 p x Pi (17) 

where Y t = y i + e; and p is the order of the equation. Equation (17) can be expressed 
more generally in matrix form as 

Y = x/3 

where Y is a (n x 1) vector of responses with error, x is a [n x (p + 1)] expanded 
model matrix, and (3 is a [( p + 1) x 1] vector of regression coefficients. In matrix 
form, OLS minimizes 

/3 = min (y - y)' (y - y) = min (V ^ x/3) ' (V - x/3^ 

The minimum is found by differentiating with respect to (3 and setting the resulting 
equation to zero or 

ll[( Y -x/j)'( Y -x/j)]=0. 

Performing the indicated differentiation yields 

— 2x ; Y + 2 (x ; x) (3 = 0. 

Therefore, the least squares normal equations are 


(x'x) (3 = x'Y. 
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Finally, solving for (3 gives 

/3 = (x'x) _1 x'Y (18) 

where x must be of full column rank. The matrix (x'x) -1 is commonly known as the 
variance- covariance matrix of the estimated regression coefficients (Myers, 1990). The 
elements of this matrix are important to the statistical properties of the estimated 
coefficients. Furthermore, the failure of x to be of full column rank leads to a matrix 
that cannot be inverted; namely, (x'x). 

4.2.2 ORTHOGONAL AND MODIFIED LEAST SQUARES 

At this point, the OLS estimator derived in the previous section addresses the 
case when the variance ratio is zero. Next, a more general estimator is derived that 
can be employed when the variance ratio is greater than zero. Based on the results 
of the simple linear model, the new estimator is expected to perform similarly. The 
approach used to derive the MLS method in the simple linear model is expanded 
upon to derive both MLS and OrthLS for the simple quadratic model. The only 
distinction between the two methods is that the variance ratio is set equal to one for 


OrthLS. 
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FIG. 14: Distance Minimized by Modified and Orthogonal Least Squares for the 
Simple Polynomial Model 


Consider the point (aq,Yi) in Figure 14. The solid curve in the figure represents 
the estimated curve, which is of the form 

1 = Po + P±X + . 

As mentioned previously, there are restrictions on how large fin can be relative to 
Pi. In normalized units, fin is restricted to values of -0.2 to 0.2. This assumption 
is typically valid in measurement system applications. The vertical line segment, a, 
represents the one-dimensional distance between the data point and the estimated 
curve, which is the distance minimized in OLS. The point ^aq , Y,* j is located where 
the line segment a intersects the estimated curve. At the point ^oq, Yf^j , the deriva- 
tive of the estimated curve is evaluated. For a simple quadratic model, the derivative 
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dy (xi) 


— A 1 + 2/3nXi, 


which is the local slope of the tangent line. This slope determines the angle, q f>, given 


= \~ tan 1 ( l,!l ii'. ) = ^ - tan 1 (fa + 2$ U X ! ) . 


It is clear that <f> in the quadratic model is a function of the location in the design 
space. From the line segment d and the slope of the tangent line, the angle a is 


defined as 


7T ( w -i ( dy(x i) 
a = 5 + (l- T )tan I 


which is the same first-order approximation used in the simple linear case. The 
distinction between MLS and OrthLS is due to the angle a being a function of the 
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variance ratio 7 . For OrthLS, 7 = 1; thus, a = —. The value of a in MLS varies 
between — and — + tan ^ 1 (fa + 2fa 1 x 1 'j. Exploiting the geometry of the problem, 
the Law of Sines yields 

sin (f) 

d = — — -a 


where a can be written as Y\ — fa — fax\ — fa\x'\. The MLS estimates are obtained 
by minimizing the squared length of d or 


Ah A) A 


E t 2 • b 

d, = nun > a, 

* ^ sin on 


For the OrthLS estimates, Equation (19) simplifies to 


n n 

A>> A, An = min X A 2 = min X t( sin A) a * 


since sin — — 1. While it is possible to solve this equation analytically, the solutions 
to Equations (19) and (20) are obtained numerically. The OLS estimates are used 
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as an initial guess for the algorithm, and since the initial guess is generally close, 
convergence of the solution occurs rather quickly. 

4.3 SIMULATION STUDY FOR THE SIMPLE QUADRATIC 
MODEL 

In the simple quadratic model, the es are assumed to be independently and iden- 
tically distributed as normal with a mean of zero and a constant variance of 
Under these conditions, the Gauss-Markov theorem holds, and the statistical prop- 
erties of the OLS estimators are well-known. Specifically, the estimated coefficients 
are unbiased, or E = (3, and the coefficients have minimal variance. It is shown 
that 

Var — a'l (x'x) -1 

where (x'x) -1 is known as the variance-covariance matrix of the /3’s. The variances 
of the estimated coefficients are located along the diagonal with the covariances 
appearing on the off-diagonal (Myers, 1990). To understand the statistical properties 
of the OrthLS and MLS estimators, a simulation was conducted to make inferences 
on the effects of MEs on the prediction capabilities of each method for the simple 
quadratic model. The prediction capabilities are again based on a variance ratio- 
weighted distance between any data point and the estimate curve. A representative 
calibration experiment is designed with n = 6 design points in order to fit the simple 
quadratic model 

V = A) + Pix + Purf + e - 

In the presence of ME, this model becomes 


y = (3 q + (3i (x + u) + fin (xi + Ui f + e. 


( 21 ) 
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The following assumptions are made about Equation (21): 

• (3q = 0, (3i = 1, and fin = 0.2. These values are representative of real mea- 
surement systems, where the sensitivity is the dominant effect. The curvature 
is limited to 20 percent of F.S. effect. 

• The errors u and e are independent and identically distributed as normal with 
means of zero and constant variances of a ^ and of, respectively. 

• <j\ and of are proportionally related through the variance ratio, 7 . 

• Wi is the error-prone value of ay. Therefore, W t and Y t are jointly distributed 
as a bivariate normal distribution 


[Wi, Yi\ ~N[( Xi ,l3 t 0 + (3 lXi + fin*!) , S] 


where E = diag (ycr/, cr/). Both Wi and Yj are independent random variables. 

The structure of the simulation is similar to the simulation executed for the simple 
linear model. From the n design points and corresponding responses, the three 
model coefficients are estimated and the MSE is calculated for each of the estimation 
methods. Five (5) confirmation points located at x = —2/3, —1/3, 0, 1/3, and 2/3 
are used to generate confirmation data to further examine the prediction capabilities 
of each estimated model. For a given calibration design, each method estimates a 
set of coefficients M times. Each set of estimated coefficients are used L times with 
the confirmation points and the MSE is calculated across (M x L) observations. 

The simulation study considered two of the three designs that were used for 
the simple linear simulation and are shown in Table 14. As mentioned previously, 
Design ^3 used in the simple linear simulation does not allow estimation of the 
curvature since there are only two unique levels of x. The features of each design are 
briefly reviewed for convenience. In Design #1, the replicates are evenly dispersed 
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TABLE 14: Designs Considered for the Simple Quadratic Simulation Study 


Design 

Design Points 

Comments 

Design #1 
Design #2 

-1,-1, 0, 0, 1,1 
-1, 0, 0, 0, 0, 1 

Replicates throughout Design Space 
Suspicion of Significant Lack-of-fit 


among the unique levels, which provides a more appropriate metric for assessing 
the constant variance assumption. Design ^2 assumes constant variance and simply 
uses the estimate of repeatability at the center (i.e. x — 0) to quote the repeatability 
throughout the design space. Table 3 shows the values of variance ratio and response 
uncertainty that are used as inputs to the simulation. The simulation employs a 
factorial combination of the variance ratio and response uncertainty. 

4.3.1 SIMULATION RESULTS 

Table 15 shows the mean estimates of the model coefficients for the two calibration 
designs and the three estimation methods. LInlike the simple linear simulation results, 
all three model coefficients are provided. From Table 15, the following inferences are 
made based on the structure of the simulation: 

• The zero-intercept estimate varies between the three estimation methods de- 
pending on the level of variance ratio and response uncertainty. In the simple 
linear simulation, all three methods yielded the same estimate of the zero- 
intercept. The impact of the method on the estimates of the coefficients is 
greater for a quadratic model than a linear model. 

• For a variance ratio of one (7 = 1), the estimates of all three model coefficients 
obtained from MLS and OrthLS are equivalent as expected. In the limit as 
7 — » 0, it is evident that the estimate of ( 3 \ from MLS converges towards the 


OLS estimate. 
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• For an uncertainty of 1 percent or better (<r e < 0.01), the three methods yield 
numerically-similar estimates of the sensitivity coefficient to four significant 
figures. 

• Larger differences between the estimates are observed when the uncertainty 
is larger than 5 percent ( cr e > 0.05). The magnitude of the difference varies 
between designs and levels of uncertainty. 

• For any level of variance ratio or uncertainty, 





OrthLS 


However, the table reveals that a similar conclusion cannot be made for the 
other two coefficients. 

Table 16 presents the variance of the mean estimates of the model coefficients for the 
two designs and three methods. Based on the results, it is clear that the variances 
are a function of the calibration design. In Design #1, the replicated points are 
evenly dispersed throughout the design space, which benefits the sensitivity estimate. 
Therefore, a more precise (less variable) estimate is obtained. Design ^2 is more 
centrally weighted in its replication strategy, which favors the zero-intercept estimate. 
Comparing the results for Design #2, it is evident that the design yields a more 
precise estimate of the zero-intercept. Between-method variability is small and from a 
practical standpoint, the variability in the MLS estimates is approximately equivalent 
to the variability in the OLS estimates. 



TABLE 15: Mean Estimates of the Model Coefficients for the Simple Quadratic Simulation 


68 



CO 

kO 

kO 


b- 

CO 

00 

CO 


00 

CO 

CO 

00 

CO 


CM 

kO 

o 

t-H 

05 

05 

05 

CM 

co 

b- 

CO 


h-3 

CO 

00 

T— 1 

co 

CO 

o 


00 

00 

CM 

CM 

CO 


05 

co 


t-H 

00 

CO 

CO 

o 

00 

t-H 

b- 

o 


r^H 

o 


05 

T — 1 


CM 

T— 1 

CM 

00 

CM 

CO 

o 

o 

CO 

t-H 

t-H 

t-H 

CM 

o 

o 

00 

o 

t-H 

05 

o 


4-3 

o 

CD 

T— 1 

o 

CD 

CM 

o 

o 

t-H 

O 

CD 

CM 

o 

o 

CM 

O 

CD 

CM 

CD 

CD 

t-H 

o 

CD 

t-H 

o 


o 

o 

T— 1 

CD 

CD 

T— 1 

CD 

o 

t-H 

CD 

o 

t-H 

CD 

o 

t-H 

CD 

CD 

t-H 

CD 

CD 

t-H 

CD 

CD 

t-H 

CD 

CD 

CM 
























































kO 

kO 


CO 

b- 

CO 

CO 

CO 

CM 

00 

b- 

00 

05 

co 

05 

kO 

o 

t-H 

05 


b- 

CM 

t-H 

kO 

CO 

3 

CO 

CO 

00 

T— 1 

CO 

o 

b- 

LO 

CO 

t-H 

CM 

CO 

kO 


o 

kO 


t-H 

00 

CO 

CM 

00 

00 



o 

m bD 


o 


05 

T— 1 

CM 

o 

t-H 

05 

b- 

CM 

o 

05 

o 

o 

05 

t-H 

t-H 

CM 

o 

O 

b- 

o 

o 

05 

o 



o 

o 

T— 1 

o 

o 

CM 

o 

05 

T— 1 

o 

o 

t-H 

o 

o 

t-H 

O 

o 

CM 

o 

o 

t-H 

o 

o 

t-H 

CD 

CD 

Q 


o 

T— 1 

CD 

CD 

T— 1 

CD 

o 

CD 

CD 

o 

t-H 

o 

o 

t-H 

CD 

CD 

t-H 

CD 

CD 

t-H 

CD 

CD 

t-H 

CD 

CD 



00 

05 

00 

CO 

05 

O 

CO 

b- 

05 

00 

t-H 

00 

05 

CO 

05 


CM 

b- 

05 


O 

CM 

05 


CO 


CO 

CO 


00 

CO 

CO 


kO 

00 

O 

CM 

CO 

kO 


o 

kO 


CO 


CO 

00 

00 

00 

CM 


o 



o 

t— 1 

b- 

T— 1 

o 

o 

t-H 

00 

b- 

CM 

o 

05 

o 

o 

05 

t-H 

o 

CM 

o 

05 

b- 

o 

O 

05 

o 


o 

o 

o 

T— 1 

o 

o 

CM 

o 

05 

t-H 

o 

CD 

t-H 

o 

o 

t-H 

O 

o 

CM 

o 

05 

t-H 

o 

CD 

t-H 

o 



o 

T— 1 

CD 

CD 

T— 1 

CD 

o 

CD 

CD 

o 

t-H 

CD 

o 

t-H 

CD 

CD 

t-H 

CD 

CD 

CD 

CD 

CD 

t-H 

CD 

CD 


CO 

00 

CO 

O 

o 

CO 

00 

CO 

t-H 

t-H 

00 

CM 

CM 

CO 

CM 

00 

05 

kO 

b- 


O 

CM 


00 

O 




CO 

00 


o 

05 

T— 1 

t-H 

kO 

00 

b- 

t-H 

CM 

co 

t-H 

00 

05 

t-H 

kO 

o 



kO 

CO 

O 



03 

o 

CM 

H 

CM 

o 

05 

t-H 

o 

05 

CO 

t-H 

CO 

t-H 

CM 

t-H 

O 

t-H 

00 

o 

o 

O 

o 

o 

O 

o 


4-3 

j_l 

o 

o 

CM 

o 

CD 

T— 1 

o 

o 

t-H 

o 

o 

t-H 

o 

O 

CM 

o 

o 

t-H 

o 

o 

CM 

o 

o 

CM 

o 


o 

o 

T— 1 

CD 

CD 

T— 1 

CD 

CD 

t-H 

CD 

CD 

t-H 

CD 

CD 

t-H 

CD 

CD 

t-H 

CD 

CD 

t-H 

CD 

CD 

t-H 

CD 

CD 

T— 1 





























00 

CO 

O 

b- 

CO 

00 

kO 

00 

05 

05 

CM 

b- 

CO 

00 

kO 

05 

kO 

b- 

CM 


CM 

00 

o 

t-H 

CO 


CO 

CO 

00 


CM 

o 

00 

05 

t-H 

kO 

CM 

CO 

kO 

o 

CO 

co 

05 

t-H 

kO 

CM 

CM 

b- 

00 

t-H 


t-H 

hn 


o 

CM 

r— 1 

T — 1 

o 

05 

o 

05 

05 

CO 

05 

CO 

CM 

o 

CM 

o 

t-H 

00 

O 

O 

O 

o 

o 

o 

o 

OAJ 


o 

o 

CM 

o 

CD 

T— 1 

o 

05 

t-H 

o 

05 

t-H 

o 

CD 

CM 

o 

CD 

t-H 

o 

o 

CM 

o 

o 

CM 

o 

c n 
CD 


o 

T— 1 

CD 

CD 

T— 1 

CD 

o 

CD 

CD 

CD 

CD 

CD 

CD 

t-H 

CD 

CD 

t-H 

CD 

CD 

t-H 

CD 

CD 

t-H 

CD 

CD 

p 





























CO 

CM 

o 

CO 

05 

t-H 


b- 

05 

t-H 

00 

CO 

CO 

00 

kO 

o 

o 


kO 

CO 

CO 

00 


t-H 

CO 


CO 

'nF 

CM 

o 


CM 

kO 

t-H 

05 

CM 

CO 

kO 

kO 

o 

CO 

CO 

co 

00 

O 

CM 

o 

b- 

00 

o 


t-H 


i_3 

T— 1 

T— 1 

CM 

T— 1 

05 

05 

t-H 

00 

05 

co 

05 

CO 

CM 

o 

CM 

o 

o 

05 

O 

o 

o 

o 

o 

o 

o 


o 

o 

CD 

CM 

CD 

05 

T— 1 

CD 

05 

t-H 

o 

05 

t-H 

o 

CD 

CM 

o 

o 

t-H 

o 

CD 

CM 

o 

CD 

CM 

o 



o 

T— 1 

CD 

o 

CD 

CD 

CD 

o 

CD 

o 

CD 

CD 

CD 

t-H 

CD 

o 

t-H 

CD 

CD 

t-H 

CD 

CD 

t-H 

CD 

CD 


4-3 



























£ 



























CD 


























CD 



























"O 

*c 

o 

1 — 1 

1 

o 


i - 1 

o 


T 1 

o 


T 1 

o 

, | 


o 

, { 

1 

o 

T 1 


o 


1 

o 

O 

ifi 

QQ_<C!Q_ 

<0$ 

<ca 

<crx 

<cj5 

<ca 

<crx 


<ca 

<crx 



<crx 









S 

0) 

o 



























o 


























kOi 



























+3 



























.s 



























’<£ 

vu 

CM 



CM 



CM 



CM 



CM 



t-H 



t-H 



t-H 



t-H 

4-3 
?— 1 
CD 

b 

o 



CD 



O 



O 



O 



CD 



CD 



CD 



CD 

CD 

















































































CD 

CD 








kO 






t-H 









kO 





_h 




kO 



CM 



t-H 



o 






kO 



CM 



t-H 

.2 

CD 

T— | 



CM 



CO 



o 



o 



t— H 



CM 



CO 



o 


"c3 




o 



o 



o 



o 






CD 



o 



CD 

£ 

Ptf 







CD 






CD 









CD 




Run 

No 

T— 1 

CM 

CO 


kO 

CO 

b- 

00 

05 



TABLE 15 - Continued 

Run Variance Uncertainty, Model Design #1 Design #2 

No Ratio, 7 Coefficient OLS MLS OrthLS OLS MLS OrthLS 
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TABLE 16: Variance of the Mean Estimates of the Model Coefficients for the Simple Quadratic Simulation 
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Table 17 shows the average MSE across all the design and confirmation points 
for the two designs and the three methods. Based on the values shown, the l-er 
estimated accuracy, in percent, is calculated using 


„ \/Mean Squared Error 

a = x 100. 

Full-Scale Range 


For Runs 11 through 15, the estimated 1-cr accuracy ranges from 2.6 to 3.6 percent 
of F.S., which is smaller than the response error of 5 percent. While Table 17 shows 
the absolute MSE values, Table 18 gives the difference between OLS and either 
MLS or OrthLS. From this table, it is once again observed that MLS has a smaller 
MSE than OLS for every combination of variance ratio and response uncertainty 
studied. Meanwhile, OrthLS only shows improvement when the variance ratio is 
one. When the variance ratio is one, OrthLS is equivalent to MLS by definition. The 
improvement in MSE is expressed in terms of percent in Table 19. There appears 
to be additional benefit in using MLS for the simple quadratic model when MEs 
are present. However, there also appears to be a larger penalty for using OrthLS as 
the variance ratio approaches zero. An analysis of variance (ANOVA) is performed 
to test whether differences in MSE are statistically detectable. From the results of 
the ANOVA shown in Tables 20 and 21, it is inferred that there is no detectable 
difference in MSE between the three methods for this simulation. This conclusion is 
also valid when looking at the hypothesis tests simultaneously. 

Table 22 relates the improvement in MSE to the known response uncertainty, cr e . 
For cases where OLS had a smaller MSE, the percent improvement is expressed as 
zero in the table. The results in the table show that MLS has a significant impact 
on the reduction in MSE for most of the cases tested during the simulation. For 
example, a 0.14 percent improvement in MSE is achieved on a 1 percent uncertainty 
measurement system using the first calibration design and a variance ratio of 0.0625. 
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This translates to an overall reduction of 14 percent, which is considered practically 
significant. Based on the results from the simple quadratic simulation, it is recom- 
mended that MLS be employed instead of OLS when the variance ratio is greater 
than 0.0625. Below a variance ratio of 0.0625, the practical benefit of MLS is small. 



TABLE 17: Mean Squared Error for the Simple Quadratic Simulation 
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TABLE 18 : Improvement in Mean Squared Error over OLS for the Simple Quadratic 
Simulation 


Run No. 

Variance 
Ratio, 7 

Uncertainty, 

07 

Design #1 
MLS OrthLS 

Design #2 
MLS OrthLS 

1 

1 

0.2 

0.0004 

0.0004 

0.0002 

0.0002 

2 

0.25 

0.2 

0.0002 

0.0000 

0.0001 

-0.0001 

3 

0.0625 

0.2 

0.0000* 

- 0.0008 

0.0000* 

-0.0010 

4 

0.01 

0.2 

0.0000* 

- 0.0008 

0.0000* 

- 0.0006 

5 

0.0001 

0.2 

0.0000* 

-0.0012 

0.0000* 

- 0.0014 

6 

1 

0.1 

0.0000* 

0.0000* 

0.0000* 

0.0000* 

7 

0.25 

0.1 

0.0000* 

0.0000 

0.0000* 

0.0000 

8 

0.0625 

0.1 

0.0000* 

0.0000 

0.0000* 

0.0000 

9 

0.01 

0.1 

0.0000* 

-0.0001 

0.0000* 

0.0000 

10 

0.0001 

0.1 

0.0000* 

-0.0001 

0.0000* 

0.0000 

11 

1 

0.05 

0.0000* 

0.0000* 

0.0000* 

0.0000* 

12 

0.25 

0.05 

0.0000* 

0.0000 

0.0000* 

0.0000 

13 

0.0625 

0.05 

0.0000* 

0.0000 

0.0000* 

0.0000 

14 

0.01 

0.05 

0.0000* 

0.0000 

0.0000* 

0.0000 

15 

0.0001 

0.05 

0.0000* 

0.0000 

0.0000* 

0.0000 

16 

1 

0.01 

0.0000* 

0.0000* 

0.0000* 

0.0000* 

17 

0.25 

0.01 

0.0000* 

0.0000 

0.0000* 

0.0000 

18 

0.0625 

0.01 

0.0000* 

0.0000 

0.0000* 

0.0000 

19 

0.01 

0.01 

0.0000* 

0.0000 

0.0000* 

0.0000 

20 

0.0001 

0.01 

0.0000* 

0.0000 

0.0000* 

0.0000 

21 

1 

0.001 

0.0000* 

0.0000* 

0.0000* 

0.0000* 

22 

0.25 

0.001 

0.0000* 

0.0000 

0.0000* 

0.0000 

23 

0.0625 

0.001 

0.0000* 

0.0000 

0.0000* 

0.0000 

24 

0.01 

0.001 

0.0000* 

0.0000 

0.0000* 

0.0000 

25 

0.0001 

0.001 

0.0000* 

0.0000 

0.0000* 

0.0000 


Note: Asterisk represents an improvement of less than 0.00005 
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TABLE 19 : Percent Improvement in Mean Squared Error over OLS for the Simple 
Quadratic Simulation 


Run No. 

Variance 
Ratio, 7 

Uncertainty, 

Design #1 
MLS OrthLS 

Design #2 
MLS OrthLS 

1 

1 

0.2 

3 . 95 % 

3 . 95 % 

2 . 59 % 

2 . 59 % 

2 

0.25 

0.2 

1 . 35 % 

- 0 . 30 % 

0.66% 

- 0 . 83 % 

3 

0.0625 

0.2 

0 . 09 % 

- 4 . 44 % 

0 . 06 % 

- 5 . 22 % 

4 

0.01 

0.2 

0.00%* 

- 3 . 73 % 

0.00%* 

- 3 . 34 % 

5 

0.0001 

0.2 

0.00%* 

- 6 . 09 % 

0.00%* 

- 6 . 17 % 

6 

1 

0.1 

0 . 97 % 

0 . 97 % 

0 . 67 % 

0 . 67 % 

7 

0.25 

0.1 

0.20% 

- 0 . 14 % 

0 . 15 % 

- 0 . 05 % 

8 

0.0625 

0.1 

0 . 03 % 

-1.10% 

0.01% 

- 0 . 63 % 

9 

0.01 

0.1 

0.00%* 

- 1 . 47 % 

0.00%* 

- 0 . 75 % 

10 

0.0001 

0.1 

0.00%* 

- 1 . 47 % 

0.00%* 

- 0 . 78 % 

11 

1 

0.05 

0 . 31 % 

0 . 31 % 

0 . 14 % 

0 . 14 % 

12 

0.25 

0.05 

0 . 08 % 

- 0 . 05 % 

0 . 04 % 

-0.01% 

13 

0.0625 

0.05 

0.00%* 

- 0 . 17 % 

0.00%* 

- 0 . 15 % 

14 

0.01 

0.05 

0.00%* 

- 0 . 33 % 

0.00%* 

-0.22% 

15 

0.0001 

0.05 

0.00%* 

- 0 . 30 % 

0.00%* 

-0.21% 

16 

1 

0.01 

0.02% 

0.02% 

0.00%* 

0.00%* 

17 

0.25 

0.01 

0.00%* 

0.00% 

0.00%* 

0.00% 

18 

0.0625 

0.01 

0.00%* 

-0.01% 

0.00%* 

0.00% 

19 

0.01 

0.01 

0.00%* 

-0.01% 

0.00%* 

-0.01% 

20 

0.0001 

0.01 

0.00%* 

-0.01% 

0.00%* 

-0.01% 

21 

1 

0.001 

0.00%* 

0.00%* 

0.00%* 

0.00%* 

22 

0.25 

0.001 

0.00%* 

0.00% 

0.00%* 

0.00% 

23 

0.0625 

0.001 

0.00%* 

0.00% 

0.00%* 

0.00% 

24 

0.01 

0.001 

0.00%* 

0.00% 

0.00%* 

0.00% 

25 

0.0001 

0.001 

0.00%* 

0.00% 

0.00%* 

0.00% 


Note: Asterisk represents an improvement of less than 0.005% 
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TABLE 20: Analysis of Variance of the Mean Squared Error for the Simple Quadratic 
Simulation - Design #1 


Source 

DF 

SS 

MS F P 

Method (OLS, MLS, OrthLS) 

2 

0.0000014 

0.0000007 0.02 0.985 

Error 

447 

0.0209744 

0.0000469 

Total 

449 

0.0209759 



TABLE 21: Analysis of Variance of the Mean Squared Error for the Simple Quadratic 
Simulation - Design ^2 


Source 

DF 

SS 

MS 

F 

P 

Method (OLS, MLS, OrthLS) 
Error 

2 

447 

0.0000016 

0.0339253 

0.0000008 

0.0000759 

0.01 

0.990 

Total 

449 

0.0339269 
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TABLE 22 : Percent Improvement in Mean Squared Error Relative to Response 
Uncertainty for the Simple Quadratic Simulation 


Run 

Variance 

Uncertainty, 

Design #1 

Design #2 

No 

Ratio, 7 


MLS 

OrthLS 

MLS 

OrthLS 

1 

1 

0.2 

10 . 56 % 

10 . 56 % 

7 . 87 % 

7 . 87 % 

2 

0.25 

0.2 

7 . 35 % 

0.00% 

5 . 05 % 

0.00% 

3 

0.0625 

0.2 

2 . 09 % 

0.00% 

1 . 73 % 

0.00% 

4 

0.01 

0.2 

0 . 41 % 

0.00% 

0 . 23 % 

0.00% 

5 

0.0001 

0.2 

0.00%* 

0.00% 

0.00%* 

0.00% 

6 

1 

0.1 

5 . 19 % 

5 . 19 % 

3 . 92 % 

3 . 92 % 

7 

0.25 

0.1 

2 . 71 % 

0.00% 

2 . 43 % 

0.00% 

8 

0.0625 

0.1 

1 . 05 % 

0.00% 

0 . 82 % 

0.00% 

9 

0.01 

0.1 

0 . 25 % 

0.00% 

0 . 13 % 

0.00% 

10 

0.0001 

0.1 

0.00%* 

0.00% 

0.00%* 

0.00% 

11 

1 

0.05 

2 . 96 % 

2 . 96 % 

1 . 72 % 

1 . 72 % 

12 

0.25 

0.05 

1 . 75 % 

0.00% 

1 . 09 % 

0.00% 

13 

0.0625 

0.05 

0 . 47 % 

0.00% 

0 . 40 % 

0.00% 

14 

0.01 

0.05 

0.10% 

0.00% 

0 . 08 % 

0.00% 

15 

0.0001 

0.05 

0.00%* 

0.00% 

0.00%* 

0.00% 

16 

1 

0.01 

0 . 77 % 

0 . 77 % 

0 . 45 % 

0 . 45 % 

17 

0.25 

0.01 

0 . 29 % 

0.00% 

0.21% 

0.00% 

18 

0.0625 

0.01 

0 . 14 % 

0.00% 

0 . 07 % 

0.00% 

19 

0.01 

0.01 

0.02% 

0.00% 

0.02% 

0.00% 

20 

0.0001 

0.01 

0.00%* 

0.00% 

0.00%* 

0.00% 

21 

1 

0.001 

0 . 05 % 

0 . 05 % 

0 . 04 % 

0 . 04 % 

22 

0.25 

0.001 

0 . 03 % 

0.00% 

0.02% 

0.00% 

23 

0.0625 

0.001 

0.00%* 

0.00% 

0.00%* 

0.00% 

24 

0.01 

0.001 

0.00%* 

0.00% 

0.00%* 

0.00% 

25 

0.0001 

0.001 

0.00%* 

0.00% 

0.00%* 

0.00% 


Note: Asterisk represents an improvement of less than 0.005% 




81 


20% 


15% 


> 

+-■ 

C 

'<5 

+-> 

o 

o 

£ 



10% 


□ <= 0.0% 

□ <=1.5% 

□ <=3.0% 

□ <= 4.5% 
■ <= 6 . 0 % 

□ <= 7.5% 
™<=9.0% 
■ >9.0% 


0.6 0.4 

Variance Ratio 


FIG. 15: Percent Improvement in Mean Squared Error of MLS for the Simple 
Quadratic Simulation - Design #1 
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FIG. 16: Percent Improvement in Mean Squared Error of MLS for the Simple 
Quadratic Simulation - Design #2 


4.4 APPLICATION OF THE MULTI-DIMENSIONAL, 
HIGHER-ORDER RESPONSE SURFACE MODEL 

Recall that an estimated model is an approximation of the true, unknown model 
and is highlighted frequently throughout this research. The models discussed to this 
point are examples of the simplest possible forms. In the context of measurement 
systems, the true behavior is often not characterized by a simple linear or quadratic 
model. This is primarily due to the fact that more than one factor often affects the 
response of the system. Unless prior knowledge exists to suggest that one of these 
simple models are sufficient, a more robust, fiexible model is recommended. In cali- 
bration applications at NASA LaRC, the second-order Taylor-series model is used for 
modeling the functional relationship of a complex measurement system. The appeal 
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of the Taylor-series model is that it is flexible, the model coefficients are estimated 
through standard regression techniques, and the model is easily expanded to include 
any number of factors of interest. Furthermore, experiments are strategically de- 
signed to ensure certain favorable properties of the estimated model. Therefore, the 
second-order Taylor-series model is the final model considered within this research. 

At NASA LaRC, most force-balances utilize the second-order Taylor-series model 
to approximate the relationship between the applied loads and the measured output 
of the strain gauges. For a standard six-component force-balance, the model is 

6 5 6 6 

Vi do + y ^ figXqj T y ^ y ^ /3ab%a,i%bi T y ' figg^a- A £ i ■ 

a= 1 a=l b=a-\- 1 a=l 

This model contains 6 linear, 15 two-factor interaction, and 6 pure quadratic terms. 
Each of the terms in the model represents some physical or electrical characteristic 
of the force-balance. For example, the first-order terms are attributed to machin- 
ing or gauging errors. The two-factor interactions are typically associated with the 
magnitude of the deflections present while the force-balance is under load (Guarino, 
1964). 

The NTF-113A/B/C are a family of single-piece force-balances used for testing in 
the NASA LaRC’s National Transonic Facility (NTF). With a full-scale normal force 
capacity of 6500 lbf . , the NTF- 113 family of force-balances are capable of measuring 
the six aerodynamic forces and moments simultaneously. These force-balances are 
gauged in the standard force configuration, which only resolves axial force and rolling 
moment directly (AIAA, 2003). Two sets of normal force and side force gauges, one 
located in the forward cage section and the other located in the aft cage section, are 
used to compute the remaining four forces and moments. The typical accuracies for 
each component for this family of force-balances is on the order of 0.1 to 0.2 percent 
of full-scale, with axial force typically having the largest amount of uncertainty. 
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The NTF-113 family of force-balances are calibrated using manual stand systems. 
Gravity-based loads, or deadweights, are applied to the force-balance through an 
intricate system of pulleys and cables. These deadweights are not known without 
error but are traceable to a NIST standard. Since the applied loads contain some 
uncertainty and because no force-balance standard exists, there are several sources of 
ME that are introduced into the calibration. Within this research, it is assumed that 
the MEs are equal across all the factors that contain MEs. However, it is possible in 
practice for the MEs to be unequal. The effect of the MEs on these complex model 
forms is studied throughout and the results are discussed shortly. 

4.5 ESTIMATION METHODS FOR THE MULTI-DIMENSIONAL, 
HIGHER-ORDER RESPONSE SURFACE MODEL 

Because the OLS estimator derived earlier for the simple quadratic model is 
also valid for multi-dimensional response surface models, it is not discussed here. 
However, the OrthLS and MLS require small modifications to the derivations in 
order to apply them to these more complex models. 

4.5.1 ORTHOGONAL AND MODIFIED LEAST SQUARES 

Recall that in the derivation of the OrthLS and MLS methods that the Law of 
Sines was used to define the vector d that was to be minimized. The Law of Sines 
required knowledge of the angles <fi and a, which were given by 

4, = * - tan-' ( 

2 \ dx J 

and 

7T . (1 w -if dy(x i)\ 
a =- + (l- 7 )tan 

Both cj) and a contain derivatives that need to be evaluated. For the simple linear 
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model, these derivatives were constants. Conversely, the derivatives for the simple 
quadratic model were a function of the location in the design space. In both cases, 
however, the geometry was relatively simple. The derivative can be interpreted as the 
rate-of-change in two-dimensional space. The corollary to this in multi-dimensional 
space is the gradient, which provides both the magnitude and direction of the rate- 
of-change of a surface. If / (aq, . . . , Xk) is the response surface function, then the 
gradient is defined as 


V/ 


' 3 / df 

dx\ 1 1 dxk 


The magnitude of the gradient 


|V/| = 




2 


replaces the derivates in <f> and a to give 



and 


7 r 


a — + (1 — 7 ) tan 1 



+ . . . + 


K 

dxi. 


The rest of the procedure for OrthLS and MLS estimation remains the same as 
previously derived. 


4.6 SIMULATION STUDY FOR THE SIX-COMPONENT 
FORCE-BALANCE 


For the six-component force-balance simulation, only one experimental design is 
considered. The design points in Table 23 are based on a spherical central composite 
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design (CCD) (Myers et al., 2009). This design is part of a larger class of second- 
order designs which are sufficient in estimating all the coefficients in the regression 
model and are typically more efficient to execute than a three-level factorial design. 
At NASA LaRC, a variant of the CCD is used for force-balance calibrations on 
the Single- Vector System (Parker et al., 2001). A full-factorial, spherical CCD in 6 
factors contains 77 unique design points, but the simulation employs a half-fraction of 
the two-level, six-factor factorial design. With a half-fraction of the factorial design, 
the CCD contains 45 unique points. In practice, fractional factorial designs are 
used as a part of screening experiments or in cases where experimental resources are 
limited. However, there are some consequences to fractionating a design, which are 
typically discussed in terms of design resolution. When estimating the coefficients 
in the second-order Taylor-series model, the fraction must be resolution V or higher. 
The half-fraction of the two-level, six-factor factorial design is a resolution VI, which 
means that there is no confounding between any of the terms in the model. Designs 
with a resolution III or IV have confounding between linear and two-factor interaction 
terms (Montgomery, 2009). 


TABLE 23: Spherical Central Composite Design for the Six-Component Force- 
Balance Simulation Study 


Design Point 

Xi 

X 2 

X3 

X 4 

x 5 

x 6 

1 

-0.6390 

-0.6390 

-0.6390 

-0.6390 

-0.6390 

-0.6390 

2 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

3 

-0.6390 

-0.6390 

0.6390 

0.6390 

0.6390 

0.6390 

4 

0.0000 

-1.0000 

0.0000 

0.0000 

0.0000 

0.0000 

5 

0.6390 

0.6390 

0.6390 

-0.6390 

-0.6390 

0.6390 

6 

-0.6390 

0.6390 

-0.6390 

0.6390 

0.6390 

0.6390 

7 

0.6390 

0.6390 

0.6390 

0.6390 

0.6390 

0.6390 

8 

0.6390 

-0.6390 

-0.6390 

0.6390 

-0.6390 

-0.6390 

9 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

1.0000 

10 

-0.6390 

-0.6390 

0.6390 

0.6390 

-0.6390 

-0.6390 

11 

0.6390 

-0.6390 

0.6390 

-0.6390 

-0.6390 

-0.6390 

12 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 



;n roint 

X\ 

13 

1.000! 

14 

0.6391 

15 

0.6391 

16 

0.6391 

17 

0.000! 

18 

-0.639 

19 

-1.000 

20 

0.000! 

21 

-0.639 

22 

0.6391 

23 

0.6391 

24 

-0.639 

25 

-0.639 

26 

0.000! 

27 

0.000! 

28 

-0.639 

29 

0.6391 

30 

0.000! 

31 

-0.639 

32 

0.6391 

33 

0.6391 

34 

0.6391 

35 

0.6391 

36 

0.000! 

37 

-0.639 

38 

0.000! 

39 

0.000! 

40 

0.000! 

41 

0.6391 

42 

0.000! 

43 

0.000! 

44 

-0.639 

45 

0.6391 

46 

0.000! 

47 

-0.639 

48 

0.000! 

49 

0.000! 

50 

-0.639 

51 

-0.639 

52 

0.000! 

53 

-0.639 

54 

0.000! 


TABLE 23 - Continued 

X2 

^3 

x 4 

0.0000 

0.0000 

0 . 000 ! 

0.6390 

-0.6390 

-0.639 

-0.6390 

-0.6390 

-0.639 

-0.6390 

0.6390 

0.6391 

0.0000 

0.0000 

0 . 000 ! 

0.6390 

0.6390 

0.6391 

0.0000 

0.0000 

0 . 000 ! 

0.0000 

0.0000 

0 . 000 ! 

■0.6390 

-0.6390 

0.6391 

0.6390 

-0.6390 

-0.639 

0.6390 

0.6390 

-0.639 

0.6390 

0.6390 

0.6391 

0.6390 

-0.6390 

0.6391 

0.0000 

0.0000 

- 1.000 

0.0000 

0.0000 

0 . 000 ! 

0.6390 

0.6390 

-0.639 

■0.6390 

-0.6390 

-0.639 

1.0000 

0.0000 

0 . 000 ! 

■0.6390 

-0.6390 

-0.639 

0.6390 

-0.6390 

0.6391 

■0.6390 

0.6390 

0.6391 

0.6390 

-0.6390 

0.6391 

■0.6390 

0.6390 

-0.639 

0.0000 

0.0000 

0 . 000 ! 

-0.6390 

0.6390 

0.6391 

0.0000 

0.0000 

0 . 000 ! 

0.0000 

0.0000 

0 . 000 ! 

0.0000 

0.0000 

1 . 000 ! 

-0.6390 

-0.6390 

0.6391 

0.0000 

0.0000 

0 . 000 ! 

0.0000 

0.0000 

0 . 000 ! 

0.6390 

0.6390 

-0.639 

0.6390 

0.6390 

0.6391 

0.0000 

0.0000 

0 . 000 ! 

0.6390 

-0.6390 

-0.639 

0.0000 

1.0000 

0 . 000 ! 

0.0000 

0.0000 

0 . 000 ! 

■0.6390 

0.6390 

-0.639 

■0.6390 

-0.6390 

0.6391 

0.0000 

0.0000 

0 . 000 ! 

0.6390 

-0.6390 

-0.639 

0.0000 

- 1.0000 

0 . 000 ! 


x 5 

x 6 

0000 

0.000! 

6390 

0.639! 

,6390 

0.639! 

,6390 

0.639! 

0000 

-1.000 

6390 

-0.639 

0000 

0.000! 

0000 

0.000! 

,6390 

0.639! 

,6390 

-0.639 

6390 

-0.639 

,6390 

0.639! 

,6390 

-0.639 

0000 

0.000! 

0000 

0.000! 

,6390 

-0.639 

6390 

-0.639 

0000 

0.000! 

6390 

0.639! 

6390 

-0.639 

6390 

-0.639 

,6390 

0.639! 

6390 

0.639! 

0000 

0.000! 

6390 

-0.639 

0000 

0.000! 

0000 

0.000! 

0000 

0.000! 

6390 

0.639! 

0000 

0.000! 

,0000 

0.000! 

6390 

0.639! 

,6390 

-0.639 

0000 

0.000! 

6390 

-0.639 

0000 

0.000! 

0000 

0.000! 

,6390 

0.639! 

6390 

-0.639 

0000 

0.000! 

,6390 

0.639! 

0000 

0.000C 






The CCD is used M = 100 times to simulate 100 calibrations of the force-balance. 


Each calibration produces a new set of estimated coefficients for each method and 
from these coefficients, the MSE based on the design points alone is calculated. For 
a given set of estimated coefficients, each set are applied L = 1000 times to the 
confirmation points shown in Table 24. The 34 confirmation points are based on a 
computer algorithm that distributes the points evenly throughout the six-dimensional 
design space. These points are used to generate confirmation data to further examine 
the prediction capabilities of each estimated model and calculate an overall MSE. In 
total, for one given set of coefficients, 34 x 1000 = 34000 confirmation points are 
applied. This simulation can be interpreted practically as the number of data points 
collected on a force-balance over the period of one calibration cycle. 


TABLE 24: Confirmation Points for the Six-Component Force-Balance Simulation 
Study 


Confirmation Point 

X\ 

X2 

x 3 

X4. 

x 5 

x 6 

1 

-0.1380 

-0.6000 

-0.4140 

0.2640 

-0.1380 

-0.6000 

2 

-0.6000 

-0.6000 

-0.6000 

0.6000 

0.6000 

0.3660 

3 

-0.6000 

-0.6000 

0.3720 

-0.6000 

-0.6000 

-0.6000 

4 

-0.1260 

0.6000 

-0.0960 

-0.6000 

-0.2280 

-0.2700 

5 

-0.1260 

0.6000 

-0.0960 

-0.6000 

-0.2280 

-0.2700 

6 

-0.6000 

-0.6000 

0.6000 

0.6000 

0.6000 

-0.6000 

7 

-0.6000 

0.6000 

0.6000 

-0.6000 

-0.6000 

0.6000 

8 

-0.5280 

0.6000 

-0.3480 

0.2880 

-0.0060 

0.6000 

9 

0.6000 

-0.6000 

-0.6000 

-0.6000 

-0.2580 

0.6000 

10 

-0.6000 

0.6000 

-0.6000 

0.6000 

-0.6000 

-0.6000 

11 

0.6000 

0.2280 

-0.5520 

0.6000 

0.6000 

0.6000 

12 

-0.5400 

-0.4080 

0.0720 

-0.4740 

0.4620 

0.6000 

13 

-0.4500 

-0.1200 

-0.6000 

-0.2880 

-0.6000 

0.1680 

14 

-0.6000 

-0.1240 

0.1500 

0.6000 

-0.2640 

-0.0480 

15 

-0.1860 

0.1740 

0.6000 

0.0000 

-0.3060 

-0.6000 

16 

-0.6000 

-0.1240 

0.1500 

0.6000 

-0.2640 

-0.0480 

17 

0.5040 

-0.5580 

-0.0942 

-0.0780 

0.6000 

-0.0840 

18 

0.6000 

0.6000 

0.6000 

-0.6000 

0.6000 

-0.6000 

19 

0.6000 

-0.1200 

-0.6000 

-0.6000 

-0.6000 

-0.6000 

20 

0.6000 

0.5280 

-0.6000 

0.2340 

0.0809 

-0.2160 

21 

0.6000 

0.6000 

0.1620 

0.6000 

-0.2270 

-0.6000 
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TABLE 24 - Continued 


Confirmation Point 

Xi 



aq 

x 5 

Xq 

22 

0.0308 

0.2520 

-0.0715 

0.6000 

0.6000 

-0.5100 

23 

0.1800 

-0.2700 

0.6000 

-0.6000 

0.1020 

0.0360 

24 

-0.1260 

-0.6000 

0.6000 

0.2390 

-0.6000 

0.6000 

25 

0.6000 

0.2160 

0.1780 

-0.2280 

-0.6000 

0.4560 

26 

0.1800 

-0.2700 

0.6000 

-0.6000 

0.1020 

0.0360 

27 

0.2880 

0.5400 

-0.6000 

-0.6000 

0.6000 

0.6000 

28 

0.6000 

-0.6000 

-0.3240 

0.0420 

-0.6000 

-0.0180 

29 

0.5040 

-0.5580 

-0.0942 

-0.0780 

0.6000 

-0.0840 

30 

0.6000 

-0.6000 

0.6000 

0.6000 

-0.6000 

-0.5340 

31 

0.3660 

-0.1800 

-0.6000 

0.6000 

-0.6000 

0.6000 

32 

-0.1860 

0.1740 

0.6000 

0.0000 

-0.3060 

-0.6000 

33 

-0.5760 

0.6000 

0.5820 

-0.0660 

0.6000 

0.0720 

34 

0.2400 

0.6000 

0.6000 

0.6000 

-0.6000 

0.2400 

35 

-0.6000 

0.5400 

-0.6000 

-0.4660 

0.1140 

0.2500 

36 

-0.6000 

0.0420 

-0.6000 

-0.6000 

0.6000 

-0.6000 

37 

0.6000 

0.6000 

0.1260 

-0.0900 

0.3360 

0.3410 

38 

0.6000 

-0.0360 

0.6000 

0.6000 

0.4020 

0.6000 


If the assumed form of the model is 

6 5 6 6 

Vi = A) + T y ^ y ^ fiabXaiXbi T y ] (3 aa% ai + Qj 

a=l a=l 6=a+l a=l 


then in the presence of ME, it becomes 


Vi = A) + ^ + U a,i) + £ £ Pab (x at + «aj + U&J 

a=l 6=a+l 


“1“ y! fiaa (%ai H“ ^di) H - Q* 
a=l 


( 22 ) 


The following assumptions are made about Equation (22) 


• The /3s are given in Table 25, which are based on an actual set of coefficients 
for the axial force component of the NTF-113C force-balance. 


• Mi = U 2 = u 3 = «4 = u$ = uq = u. The six MEs are assumed to be equal 
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and are represented by u. The errors u and e are independently and identically 
distributed as normal with means of zero and constant variances of a 2 and of, 
respectively. 

• of and of are proportionally related through the variance ratio, 7. 

• (Ilf ., ir (ii ) are the error-prone values of [x\ ir . . . ,x Xi ). Therefore, 

(ITf , . . . , Wq J and Y i are jointly distributed as a multivariate normal distribu- 
tion 


iw u ,...,w s ,,y ( ] 

~ N 


6 


6 


Xi i ,...,X 6i ,(3 0 + Yl + Y Y ftab'Kcii'X'bi E Paa^ai } ) ^ | 

d=l d=l b=d -\- 1 d=l 


where E = diag (yof , yof, yof, yof, yof, y of, of). Both (W u , . . . , W 6i ) and Y t 
are independent random variables. 
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TABLE 25: Model Coefficients for the Six-Component Force-Balance Simulation 
Study 


First-Order Terms 

Interaction Terms 

Quadratic Terms 

A 

0.0000 

A 2 

0.0321 

Ai 

-0.1015 

A 

0.1334 

A3 

0.0934 

A 2 

0.0007 

A 

1.0000 

A4 

0.0014 

A3 

0.0387 

A 

-0.0095 

A 5 

-0.0027 

Pu 

0.0624 

A 

0.0257 

A 6 

0.0023 

A5 

-0.0025 

A 

-0.0029 

A3 

-0.0025 

Pm 

-0.0299 

A 

0.0006 

A4 

0.0000 





As 

-0.0003 





A 6 

0.0009 





A4 

0.0000 





P 35 

-0.0013 





@36 

-0.0001 





A5 

-0.0683 





A6 

-0.0105 





A6 

0.0306 




4.6.1 SIMULATION RESULTS 

Table 26 shows the mean estimates of the model coefficients for two combinations 
of variance ratio and response uncertainty: 7 = 0.25, a e = 0.1 and 7 = 0.0625, 
cr e = 0.001. In the first case, larger differences in the estimates are observed in the 
linear and quadratic terms. For example, the primary sensitivity, A> ranges from 
0.9920 for OLS to 1.0284 for OrthLS. This corresponds to a 3.6 percent difference 
over the F.S. range. Smaller effects, such as An vary as much as 2.5 percent of F.S. 
It is also interesting to note the trends in the estimates for the first case are not 
consistent. Simply put, 

AoLS — A M LS — AorthLS 

does not hold for all 28 estimates. The second case of 7 = 0.0625 and a e = 0.001 
shown in the table is representative of a NTF-113 calibration. Unlike the first case, 
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these estimates are numerically equivalent to four significant figures across the three 
methods. Table 27 shows the variance of the coefficient estimates for the three 
methods. Between the methods, the variances are similar for a given coefficient. 
Furthermore, variances within like coefficients are also similar. For example, the 
variance of the linear coefficients is approximately 0.0006. These variances are again 
attributed to the design rather than the estimation method; thus, the variability in 
the MLS estimates are equal to the OLS estimates. 
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TABLE 26: Mean Estimates of the Model Coefficients for the Six-Component Force- 
Balance Simulation 


Model Term 

7 = 

0.25, = 

0.1 

7 = 0.0625, cr e = 

0.001 

OLS 

MLS 

OrthLS 

OLS 

MLS 

OrthLS 

Po 

0.0069 

0.0071 

0.0071 

0.0000 

0.0000 

0.0000 

h 

0.1343 

0.1364 

0.1388 

1.0000 

1.0000 

1.0000 

h 

0.9920 

1.0087 

1.0284 

-0.0801 

-0.0801 

-0.0801 

h 

-0.0020 

-0.0019 

-0.0018 

0.0900 

0.0900 

0.0900 

Pa 

0.0144 

0.0148 

0.0152 

-0.0300 

-0.0300 

-0.0300 

A 

-0.0104 

-0.0108 

-0.0112 

-0.0500 

-0.0500 

-0.0500 

Ps 

-0.0037 

-0.0038 

-0.0038 

0.3500 

0.3500 

0.3500 

$12 

0.0364 

0.0380 

0.0397 

0.0500 

0.0500 

0.0500 

Pl3 

0.0918 

0.0926 

0.0926 

0.0000 

0.0000 

0.0000 

PlA 

0.0203 

0.0212 

0.0221 

-0.0100 

-0.0100 

-0.0100 

Pl 5 

0.0075 

0.0084 

0.0095 

-0.0201 

-0.0201 

-0.0201 

Pl6 

0.0052 

0.0043 

0.0029 

0.0699 

0.0699 

0.0699 

P ‘23 

-0.0029 

-0.0028 

-0.0023 

-0.0500 

-0.0500 

-0.0501 

P 24 

0.0190 

0.0184 

0.0180 

-0.0099 

-0.0099 

-0.0099 

P25 

0.0007 

0.0012 

0.0020 

0.0000 

0.0000 

0.0000 

P 26 

0.0020 

0.0021 

0.0023 

0.0300 

0.0300 

0.0300 

P34 

0.0009 

0.0003 

0.0000 

0.0000 

0.0000 

0.0000 

p 35 

0.0027 

0.0027 

0.0029 

-0.0199 

-0.0199 

-0.0199 

P 36 

0.0110 

0.0105 

0.0099 

-0.1000 

-0.1000 

-0.1000 

Pa5 

-0.0671 

-0.0685 

-0.0697 

-0.0300 

-0.0300 

-0.0300 

Pas 

-0.0119 

-0.0126 

-0.0130 

-0.4500 

-0.4500 

-0.4499 

Pm 

0.0399 

0.0399 

0.0401 

0.0200 

0.0200 

0.0200 

P 11 

-0.1151 

-0.1265 

-0.1401 

0.1200 

0.1200 

0.1199 

P22 

-0.0003 

0.0015 

0.0041 

-0.0100 

-0.0100 

-0.0100 

P33 

0.0432 

0.0473 

0.0512 

0.0399 

0.0399 

0.0400 

Paa 

0.0709 

0.0785 

0.0884 

-0.0498 

-0.0498 

-0.0498 

P55 

-0.0141 

-0.0143 

-0.0151 

0.0100 

0.0100 

0.0100 

Pm 

-0.0357 

-0.0382 

-0.0404 

-0.0401 

-0.0401 

-0.0401 
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TABLE 27: Variance of the Mean Estimates of the Model Coefficients for the Six- 
Component Force-Balance Simulation 


Model Term 

7 = 

0.25, a e 

= 0.1 

7 = 0.0625, cr e 

= 0.001 

OLS 

MLS 

OrthLS 

OLS 

MLS 

OrthLS 

A) 

0.0006 

0.0006 

0.0006 

0.0000 

0.0000 

0.0000 

h 

0.0006 

0.0007 

0.0007 

0.0000 

0.0000 

0.0000 

h 

0.0006 

0.0006 

0.0006 

0.0000 

0.0000 

0.0000 

h 

0.0005 

0.0006 

0.0006 

0.0000 

0.0000 

0.0000 

$a 

0.0006 

0.0006 

0.0006 

0.0000 

0.0000 

0.0000 

$5 

0.0005 

0.0005 

0.0005 

0.0000 

0.0000 

0.0000 

$e 

0.0007 

0.0007 

0.0007 

0.0000 

0.0000 

0.0000 

$12 

0.0015 

0.0015 

0.0016 

0.0000 

0.0000 

0.0000 

$13 

0.0018 

0.0018 

0.0019 

0.0000 

0.0000 

0.0000 

$IA 

0.0015 

0.0015 

0.0016 

0.0000 

0.0000 

0.0000 

$15 

0.0020 

0.0020 

0.0020 

0.0000 

0.0000 

0.0000 

$16 

0.0019 

0.0020 

0.0020 

0.0000 

0.0000 

0.0000 

$23 

0.0017 

0.0017 

0.0017 

0.0000 

0.0000 

0.0000 

$2A 

0.0017 

0.0017 

0.0018 

0.0000 

0.0000 

0.0000 

$25 

0.0019 

0.0019 

0.0020 

0.0000 

0.0000 

0.0000 

$26 

0.0018 

0.0018 

0.0018 

0.0000 

0.0000 

0.0000 

$3A 

0.0013 

0.0013 

0.0013 

0.0000 

0.0000 

0.0000 

$35 

0.0019 

0.0019 

0.0020 

0.0000 

0.0000 

0.0000 

$36 

0.0018 

0.0019 

0.0019 

0.0000 

0.0000 

0.0000 

$A5 

0.0019 

0.0019 

0.0019 

0.0000 

0.0000 

0.0000 

$A6 

0.0019 

0.0019 

0.0019 

0.0000 

0.0000 

0.0000 

$56 

0.0016 

0.0016 

0.0016 

0.0000 

0.0000 

0.0000 

$11 

0.0049 

0.0050 

0.0052 

0.0000 

0.0000 

0.0000 

$22 

0.0051 

0.0054 

0.0063 

0.0000 

0.0000 

0.0000 

$33 

0.0043 

0.0045 

0.0048 

0.0000 

0.0000 

0.0000 

$AA 

0.0035 

0.0036 

0.0038 

0.0000 

0.0000 

0.0000 

$55 

0.0039 

0.0041 

0.0042 

0.0000 

0.0000 

0.0000 

$66 

0.0049 

0.0052 

0.0055 

0.0000 

0.0000 

0.0000 
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Table 28 shows the average MSE across all the design and confirmation points for 
each of the three methods. The case of 7 = 0.0625 and a e = 0.001 is closest to the 
actual calibration of a NTF-113 force-balance. Based on the results from this case, 
the quote accuracy is 


/4.5190 x 10- 


Full-Scale Range 


or 0.07 percent of F.S. Table 29 shows the difference in MSE between OLS and 
either MLS or OrthLS. While the MSE for MLS is smaller in all the cases studied, 
the absolute difference is once again smaller than 5.0 x 10~ 5 . In terms of percent 
difference, as shown in Table 30, there are more noticeable improvements in MSE 
for variance ratios greater than 0.25 and uncertainties larger than 0.1. From the 
ANOVA in Table 31, there still remains strong evidence ( p > 0.05) that the three 
methods do not provide a detectable difference in MSE. 

Lastly, Table 32 reiterates the practical benefits of the improvement in MSE at- 
tained by the method of MLS. For instances where the OLS had a smaller MSE, the 
percent improvement is expressed as zero. The improvement in MSE by using MLS 
estimation is once again observed in the results. For example, 0.2 percent improve- 
ment in MSE is achieved on a measurement system with a 0.1 percent uncertainty 
in the response with a variance ratio of 0.025. This is a substantial difference over 
the improvement in the same combination of uncertainty and variance ratio in the 
simple linear study, where the percent improvement was 0.02 percent. Therefore, in 
complex model forms, such as the Taylor-series model, the impact of MLS is much 
greater. It is once again concluded that the largest benefit of MLS is in situations 
where the variance ratio is greater than 0.0625. However, the data do suggest that 
the practical benefit of MLS at lower variance ratios may become more pronounced 
in higher-order models. 
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TABLE 28: Mean Squared Error over OLS for the Six-Component Force-Balance 
Simulation 


Run No. 

Variance 
Ratio, 7 

Uncertainty, a e 

OLS 

MLS 

OrthLS 

1 

1 

0.2 

0.0096 

0.0090 

0.0090 

2 

0.25 

0.2 

0.0136 

0.0134 

0.0137 

3 

0.0625 

0.2 

0.0169 

0.0169 

0.0179 

4 

0.01 

0.2 

0.0185 

0.0185 

0.0198 

5 

0.0001 

0.2 

0.0202 

0.0202 

0.0219 

6 

1 

0.1 

0.0024 

0.0024 

0.0024 

7 

0.25 

0.1 

0.0034 

0.0034 

0.0034 

8 

0.0625 

0.1 

0.0045 

0.0045 

0.0045 

9 

0.01 

0.1 

0.0048 

0.0048 

0.0048 

10 

0.0001 

0.1 

0.0050 

0.0050 

0.0051 

11 

1 

0.05 

0.0006 

0.0006 

0.0006 

12 

0.25 

0.05 

0.0008 

0.0008 

0.0008 

13 

0.0625 

0.05 

0.0011 

0.0011 

0.0011 

14 

0.01 

0.05 

0.0012 

0.0012 

0.0012 

15 

0.0001 

0.05 

0.0012 

0.0012 

0.0012 

16 

1 

0.01 

0.0000 

0.0000 

0.0000 

17 

0.25 

0.01 

0.0000 

0.0000 

0.0000 

18 

0.0625 

0.01 

0.0000 

0.0000 

0.0000 

19 

0.01 

0.01 

0.0000 

0.0000 

0.0000 

20 

0.0001 

0.01 

0.0000 

0.0000 

0.0000 

21 

1 

0.001 

0.0000 

0.0000 

0.0000 

22 

0.25 

0.001 

0.0000 

0.0000 

0.0000 

23 

0.0625 

0.001 

0.0000 

0.0000 

0.0000 

24 

0.01 

0.001 

0.0000 

0.0000 

0.0000 

25 

0.0001 

0.001 

0.0000 

0.0000 

0.0000 
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TABLE 29 : Improvement in Mean Squared Error over OLS for the Six-Component 
Force-Balance Simulation 


Run No. 

Variance Ratio, 7 

Uncertainty, cr e 

MLS 

OrthLS 

1 

1 

0.2 

0.0005 

0.0005 

2 

0.25 

0.2 

0.0002 

-0.0001 

3 

0.0625 

0.2 

0.0000* 

-0.0010 

4 

0.01 

0.2 

0.0000* 

- 0.0013 

5 

0.0001 

0.2 

0.0000* 

- 0.0017 

6 

1 

0.1 

0.0000* 

0.0000* 

7 

0.25 

0.1 

0.0000* 

0.0000 

8 

0.0625 

0.1 

0.0000* 

0.0000 

9 

0.01 

0.1 

0.0000* 

-0.0001 

10 

0.0001 

0.1 

0.0000* 

-0.0001 

11 

1 

0.05 

0.0000* 

0.0000* 

12 

0.25 

0.05 

0.0000* 

0.0000 

13 

0.0625 

0.05 

0.0000* 

0.0000 

14 

0.01 

0.05 

0.0000* 

0.0000 

15 

0.0001 

0.05 

0.0000* 

0.0000 

16 

1 

0.01 

0.0000* 

0.0000* 

17 

0.25 

0.01 

0.0000* 

0.0000* 

18 

0.0625 

0.01 

0.0000* 

0.0000 

19 

0.01 

0.01 

0.0000* 

0.0000 

20 

0.0001 

0.01 

0.0000* 

0.0000 

21 

1 

0.001 

0.0000* 

0.0000* 

22 

0.25 

0.001 

0.0000* 

0.0000* 

23 

0.0625 

0.001 

0.0000* 

0.0000 

24 

0.01 

0.001 

0.0000* 

0.0000 

25 

0.0001 

0.001 

0.0000* 

0.0000 


Note: Asterisk represents an improvement of less than 0.00005 
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TABLE 30 : Percent Improvement in Mean Squared Error over OLS for the Six- 
Component Force-Balance Simulation 


Run No. 

Variance Ratio, 7 

Uncertainty, cr e 

Design #1 
MLS OrthLS 

1 

1 

0.2 

5 . 48 % 

5 . 48 % 

2 

0.25 

0.2 

1.21% 

- 0 . 58 % 

3 

0.0625 

0.2 

0.11% 

- 6 . 13 % 

4 

0.01 

0.2 

0.00%* 

- 7 . 01 % 

5 

0.0001 

0.2 

0.00%* 

- 8 . 46 % 

6 

1 

0.1 

1 . 09 % 

1 . 09 % 

7 

0.25 

0.1 

0 . 27 % 

- 0 . 07 % 

8 

0.0625 

0.1 

0.02% 

- 0 . 81 % 

9 

0.01 

0.1 

0.00%* 

- 1 . 08 % 

10 

0.0001 

0.1 

0.00%* 

- 1 . 28 % 

11 

1 

0.05 

0 . 25 % 

0 . 25 % 

12 

0.25 

0.05 

0 . 07 % 

0.00% 

13 

0.0625 

0.05 

0.00%* 

-0.21% 

14 

0.01 

0.05 

0.00%* 

- 0 . 27 % 

15 

0.0001 

0.05 

0.00%* 

- 0 . 27 % 

16 

1 

0.01 

0.01% 

0.01% 

17 

0.25 

0.01 

0.00%* 

0.00%* 

18 

0.0625 

0.01 

0.00%* 

-0.01% 

19 

0.01 

0.01 

0.00%* 

-0.01% 

20 

0.0001 

0.01 

0.00%* 

-0.01% 

21 

1 

0.001 

0.00%* 

0.00%* 

22 

0.25 

0.001 

0.00%* 

0.00%* 

23 

0.0625 

0.001 

0.00%* 

0.00% 

24 

0.01 

0.001 

0.00%* 

0.00% 

25 

0.0001 

0.001 

0.00%* 

0.00% 


Note: Asterisk represents an improvement of less than 0.005% 


TABLE 31 : Analysis of Variance of the Mean Squared Error for the Six-Component 
Force-Balance Simulation 


Source 

DF 

SS 

MS 

F 

P 

Method (OLS, MLS, OrthLS) 
Error 

2 

4047 

0.0000253 

0.2165413 

0.0000126 

0.0000535 

0.24 

0.790 

Total 

4049 

0.2165666 
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TABLE 32 : Percent Improvement in Mean Squared Error Relative to Response 
Uncertainty for the Six-Component Force-Balance Simulation 


Run No. 

Variance Ratio, 7 

Uncertainty, a e 

Design #1 
MLS OrthLS 

1 

1 

0.2 

11 . 45 % 

11 . 45 % 

2 

0.25 

0.2 

6 . 41 % 

0.00% 

3 

0.0625 

0.2 

2 . 16 % 

0.00% 

4 

0.01 

0.2 

0 . 37 % 

0.00% 

5 

0.0001 

0.2 

0.00%* 

0.00% 

6 

1 

0.1 

5 . 16 % 

5 . 16 % 

7 

0.25 

0.1 

3 . 06 % 

0.00% 

8 

0.0625 

0.1 

1 . 04 % 

0.00% 

9 

0.01 

0.1 

0 . 18 % 

0.00% 

10 

0.0001 

0.1 

0.00%* 

0.00% 

11 

1 

0.05 

2 . 46 % 

2 . 46 % 

12 

0.25 

0.05 

1 . 57 % 

0.00% 

13 

0.0625 

0.05 

0 . 52 % 

0.00% 

14 

0.01 

0.05 

0 . 09 % 

0.00% 

15 

0.0001 

0.05 

0.00%* 

0.00% 

16 

1 

0.01 

0 . 60 % 

0 . 60 % 

17 

0.25 

0.01 

0 . 35 % 

0 . 05 % 

18 

0.0625 

0.01 

0.12% 

0.00% 

19 

0.01 

0.01 

0.02% 

0.00% 

20 

0.0001 

0.01 

0.00%* 

0.00% 

21 

1 

0.001 

0 . 32 % 

0 . 32 % 

22 

0.25 

0.001 

0.20% 

0.02% 

23 

0.0625 

0.001 

0 . 07 % 

0.00% 

24 

0.01 

0.001 

0.01% 

0.00% 

25 

0.0001 

0.001 

0.00%* 

0.00% 


Note: Asterisk represents an improvement of less than 0.005% 
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FIG. 17: Percent Improvement in Mean Squared Error of MLS for the Six-Component 
Force-Balance Simulation 
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CHAPTER 5 


CONCLUSIONS 


Measurement errors (ME) are a source of variability in a calibration experiment that 
are often ignored due to the complexities associated with the regression analysis. 
Standard regression techniques, such as ordinary least squares (OLS), are not appro- 
priate in situations that contain MEs since these methods assume that the factors are 
known without error. In practice, it is more appropriate to assume that the factors 
are known within a certain level of uncertainty. Methods, such as those proposed by 
Adcock (1878), Kummel (1879), and Doming (1931, 1964), are available for instances 
where MEs are present in the data, but these methods are under-developed or limited 
in calibration applications at NASA. For example, there are no known methods that 
are capable of estimating the model for a 6-component force-balance while account- 
ing for any MEs. Furthermore, it is reasonable to assume that a variance ratio of 1 
is unlikely in calibration, which is the basis of using a method like orthogonal least 
squares (OrthLS). Therefore, the method of modified least squares (MLS) derived 
within this research is a method that accounts for any MEs and is employable for 
simple or complex mathematical models. 

Based on the scope of the simulations, the following are contributions of the MLS 
method: 

• New definition of residual error. In OLS, this is defined as the difference in the 
observed and predicted values of the response, or Y t — Y t . However, the residual 
error should include other sources of variability if present, such as MEs. For 
MLS, the residual error is defined as a variance ratio-weighted distance between 
the observed and predicted data point. Geometrically, this is represented by 
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the length d, which is the distance that is minimized in MLS. In terms of quoted 
accuracy of a measurement system, the former definition of residual error is an 
underestimate of the true accuracy. 

• Smaller mean squared error (MSE). The MSE for MLS is smaller than the 
MSE for OLS for all the cases studied, and depending on the variance ratio 
and accuracy of the measurement system, the improvement in MSE ranges from 
0 to 5.5 percent. The improvements are attributed to small differences in the 
estimates of the model coefficients. However, the practical impact of MLS is 
seen when the MSE is compared with the accuracy of the measurement system. 
Even in cases where the accuracy is high (i.e. a e < 0.01), the improvements in 
MSE are meaningful. 

The pressure transducers for EFT-1 are an example of a simple linear instrument 
with a stated accuracy of 0.1 percent of full-scale (F.S.). With the NIST-traceable 
uncertainty in the reference standard, the variance ratio for this type of calibration is 
on the order of 0.0625. From the results of the simulation, the practical improvement 
in MSE is small, approximately 0.01 percent. The accuracies of the six components 
on the NTF-113 family of force-balances range from 0.1 percent to 0.5 percent of 
F.S. Considering the same value of variance ratio of 0.0625, the results indicate a 
more substantial improvement in MSE by employing MLS. The improvement for the 
6-component is at least 0.07 percent. 

For multiple-component calibrations, the method of MLS is derived assuming that 
the MEs are equal across the k factors. In practice, it is unlikely that the variances 
are equal. Unequal errors translate into multiple variance ratios instead of a single 
variance ratio, as used in the 6-component force-balance simulation. By allowing 
multiple variance ratios, the individual factors are more appropriately weighted based 
on the respective uncertainties. As a result, the potential impact of MLS could be 
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greater when considering individual variance ratios. 

Model reduction is a technique that is frequently used in calibration applications 
to remove terms from the estimated model that do not detectably influence the 
response. A statistical hypothesis test is used to evaluate whether a term in the 
model is influential to a specified significance level. The test is based on Student’s 
t- distribution (Casella and Berger, 2002). The model term is considered statistically 
significant if 


Pi 


s.e. p t 


> in- 


n-2, 0.025 


where s.e. 



is the standard error if the Pi estimate, which is given by 


s.e. (p^j 



5 


and t n - 2 , 0.025 is the f-statistic with n — 2 degrees of freedom and a 95 percent con- 
fidence level. While model reduction is not discussed within the research, it is sus- 
pected that the terms deemed significant in OLS would also be significant in MLS 
and OrthLS. The same is suspected of the terms that are determined not to be sig- 
nificant. This rationale is based on the results that for any given estimated model 
coefficient, the variance of the estimate is nearly equivalent across the three methods. 

Finally, the simulations conducted within this research only represent a small 
sample of all possible simulations that are possible. The simulations are representa- 
tive of typical calibrations and life-cycles of measurement systems at NASA. However, 
to understand the general statistical properties of the MLS estimators, one of the 
following need to be completed: 

• Analytical derivation of the properties. The estimators are ratios of random 
variables, which complicates the derivation process. The delta method can be 
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used to obtain an asymptotic approximation of the ratios (Parker et ah, 2010). 

• Larger-scale simulation. In order to make general inferences based on numerical 
results, the simulation should consider a large number of additional calibration 
designs, variance ratios, and measurement system accuracies. Furthermore, the 
number of times that the coefficients are estimated in the simulation should be 
expanded since 100 is considered a small number in estimating variance. 

From a practical standpoint, a larger simulation is the logical choice in proceeding 
with understanding the statistical properties of the estimators. 
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APPENDIX A 

SIMPLE LINEAR SIMULATION CODE 


function [dist, ccdist, bmean, bvar] = simplinsim5(n,N,std, gamma) 

7. n: number of calibrations 

7 . N: number of uses per calibration 

% std: measurement system accuracy 

7. gamma: variance ratio (delta-to-epsilon) 

warning off all; 

options = optimsetp Algorithm’ , ’ levenberg-marquardt J , ... 

’ TolFun’ , le-12, ’ TolX’ , le-12, ’ Display J , ’off ’) ; 


tic ; 

7 . 7 . Design Points 

xl = [-1; -1; 0; 0; 1; 1]; 
x2 = [-1; 0; 0; 0; 0; 1] ; 
x3 = [-1; -1; -1; 1; 1; 1] ; 

xlm = [1-1; 1-1; 10; 10; 11; 11]; 
x2m = [1 -1; 1 0; 1 0; 1 0; 1 0; 1 1] ; 
x3m = [1 -1; 1 -1; 1 -1; 11; 11; 11]; 

7 . 7 . Confirmation Points 

xconf = [-1; -0.5; 0; 0.5; 1]; 

7 . 7 . True Model and Responses 
7. b0 = 0, bl = 1 


yl = xl; 
y2 = x2 ; 
y3 = x3 ; 


yconf = xconf; 
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°/„°/o Error in x and y Sampled from Bivariate Normal Distribution 

sigma = [gamma* (std~2) 0; 0 (std~2)] ; 
for i = l:n 

for j = 1:6 

cal(i) ,desl(j , : ) = mvnrnd( [xl (j , 1) yl (j , 1)] , sigma) ; 
yobs_desl ( j , i) = cal(i) .desl(j ,2) ; 

cal(i) ,des2(j , : ) = mvnrnd( [x2(j , 1) y2(j , 1)] , sigma) ; 
yobs_des2(j , i) = cal(i) ,des2(j ,2) ; 

cal(i) ,des3(j , : ) = mvnrnd( [x3(j , 1) y3(j , 1)] , sigma) ; 
yobs_des3(j , i) = cal(i) .des3(j ,2) ; 

end 

°/:/ 0 Ordinary Least Squares Estimation 

bordl(:,i) = regress (yobs_desl (:, i) ,xlm) ; 
bord2(:,i) = regress (yobs_des2( i) ,x2m) ; 
bord3(:,i) = regress (yobs_des3( i) ,x3m) ; 

°/:/ 0 Ordinary Least Squares Distance 

for j = 1:6 

phi_ordl ( j , i) = (pi/2) - atan(bordl (2, i) ) ; 
phi_ord2 ( j , i) = (pi/2) - atan(bord2(2,i)) ; 
phi_ord3(j , i) = (pi/2) - atan(bord3(2, i) ) ; 

alpha_ordl ( j , i) = (pi/2) + (1-gamma) *atan(bordl (2 , i) ) ; 
alpha_ord2(j , i) = (pi/2) + ( 1-gamma) *atan(bord2 (2, i) ) ; 
alpha_ord3(j , i) = (pi/2) + ( 1-gamma) *atan(bord3 (2, i) ) ; 

dist_ordl ( j , i) = ( (sin(phi_ordl ( j , i) ) / ... 
sin(alpha_ordl ( j , i) ) ) ... 

(yobs_desl ( j , i) - bordl(l,i) - bordl (2 , i) *xl ( j , 1) ) ) ; 
dist_ord2(j , i) = ( (sin(phi_ord2(j , i) )/ ... 
sin(alpha_ord2(j ,i)))* ... 

(yobs_des2 ( j , i) - bord2(l,i) - bord2(2,i)*x2(j , 1))) ; 
dist_ord3(j , i) = ( (sin(phi_ord3(j , i) )/ ... 
sin(alpha_ord3(j ,i)))* ... 

(yobs_des3(j , i) - bord3(l,i) - bord3(2 , i) *x3( j , 1) ) ) ; 


end 



U Orthogonal Least Squares Estimation 

borthl(:,i) = lsqnonlin(@(borthl) (sin( (pi/2) - ... 

atan(borthl (2) ) )/sin(pi/2) ) * (yobs_desl (1 : 6 , i) - ... 
borthl(l) - borthl (2) *xl (1 : 6 , 1) ) ,bordl (:, i) ,[],[], options) 

borth2(:,i) = lsqnonlin(@(borth2) (sin( (pi/2) - ... 

atan(borth2(2) ) )/sin(pi/2) ) * (yobs_des2(l : 6 , i) - ... 
borth2(l) - borth2(2)*x2(l : 6, 1) ) ,bord2( i) ,[],[], options) 

borth3(:,i) = lsqnonlin(@(borth3) (sin( (pi/2) - ... 
atan(borth3(2) ) )/sin(pi/2) ) * (yobs_des3(l : 6 , i) - ... 
borth3(l) - borth3(2) *x3(l : 6 , 1) ) ,bord3( i) ,[],[], options) 

for j = 1:6 

phi_orthl ( j , i) = (pi/2) - atan(borthl (2 , i) ) ; 
phi_orth2(j , i) = (pi/2) - atan(borth2(2,i)) ; 
phi_orth3(j , i) = (pi/2) - atan(borth3(2,i)) ; 

alpha_orthl ( j , i) = (pi/2) + (1-gamma) *atan(borthl (2 , i) ) ; 
alpha_orth2 ( j , i) = (pi/2) + (1-gamma) *atan(borth2 (2 , i) ) ; 
alpha_orth3(j , i) = (pi/2) + (l-gamma)*atan(borth3(2,i)) ; 

dist_orthl (j , i) = ( (sin(phi_orthl ( j , i) ) / ... 
sin(alpha_orthl (j , i) ) ) * (yobs_desl ( j , i) - ... 
borthl(l,i) - borthl (2 , i) *xl (j , 1) )) ; 
dist_orth2(j , i) = ( (sin(phi_orth2( j , i) ) / ... 
sin(alpha_orth2(j ,i)))*(yobs_des2(j ,i) ~ ••• 
borth2(l,i) ~ borth2(2,i)*x2(j , 1))) ; 
dist_orth3(j , i) = ( (sin(phi_orth3( j , i) ) / ... 
sin(alpha_orth3(j ,i)))*(yobs_des3(j ,i) ~ ••• 
borth3(l,i) ~ borth3(2,i)*x3(j , 1))) ; 

end 

U Modified Least Squares 

bmodl(:,i) = lsqnonlin(@(bmodl) (sin( (pi/2) - atan(bmodl(2)))/ 
sin((pi/2) + (1-gamma) *atan(bmodl (2) ))) * ... 

(yobs_desl (1 : 6 , i) ~ bmodl(l) - bmodl (2) *xl (1 : 6 , 1) ) , ... 
bordl (:, i) >[],[] > options) ; 

bmod2(:,i) = lsqnonlin(@(bmod2) (sin( (pi/2) - atan(bmod2(2)))/ 
sin((pi/2) + (1-gamma) *atan(bmod2 (2) ))) * ... 

(yobs_des2(l : 6 , i) - bmod2(l) - bmod2 (2) *x2 (1 : 6 , 1) ) , ... 
bord2 ( : , i) ,[],[] , options) ; 

bmod3(:,i) = lsqnonlin(@(bmod3) (sin( (pi/2) - atan(bmod3(2)))/ 
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sin((pi/2) + (l-gairnna)*atan(bmod3(2))))* ... 
(yobs_des3(l : 6 , i) - bmod3(l) - bmod3(2) *x3(l : 6 , 1) ) , ... 
bord3 (:, i) ,[],[] , options) ; 

for j = 1:6 

phi_raodl ( j , i) = (pi/2) - atan(bmodl (2, i) ) ; 
phi_raod2 ( j , i) = (pi/2) - atan(bmod2(2,i)) ; 
phi_raod3(j , i) = (pi/2) - atan(bmod3(2, i) ) ; 

alpha_modl ( j , i) = (pi/2) + (1-gamma) *atan(bmodl (2 , i) ) ; 
alpha_mod2(j , i) = (pi/2) + ( 1-gamma) *atan(bmod2 (2, i) ) ; 
alpha_mod3(j , i) = (pi/2) + ( 1-gamma) *atan(bmod3 (2, i) ) ; 

dist_modl ( j , i) = ( (sin(phi_modl ( j , i) ) / ... 

sin(alpha_modl ( j , i) ) ) * (yobs_desl ( j , i) - ... 
bmodl(l,i) - bmodl (2, i) *xl (j , 1) ) ) ; 
dist_mod2(j , i) = ( (sin(phi_mod2( j , i) ) / ... 

sin(alpha_mod2(j , i) ) ) * (yobs_des2 ( j ,i) - ... 
bmod2(l,i) - bmod2 (2, i) *x2( j , 1) ) ) ; 
dist_mod3(j , i) = ( (sin(phi_mod3(j , i) )/ ... 

sin(alpha_mod3(j ,i)))*(yobs_des3(j ,i) - ... 
bmod3(l,i) - bmod3(2, i) *x3(j , 1) ) ) ; 

end 

end 

for j = 1:6 

dist_ordl ( j ,n+l) = 0; 
dist_ord2 ( j ,n+l) = 0; 
dist_ord3(j ,n+l) = 0; 

dist_modl ( j ,n+l) = 0; 
dist_mod2 ( j ,n+l) = 0; 
dist_mod3(j ,n+l) = 0; 

dist_orthl ( j ,n+l) = 0; 
dist_orth2(j ,n+l) = 0; 
dist_orth3(j ,n+l) = 0; 

dist_ordl ( j ,n+2) = mean(dist_ordl (j , 1 :n) ) ; 
dist_ord2 ( j ,n+2) = mean(dist_ord2(j , 1 :n) ) ; 
dist_ord3(j ,n+2) = mean(dist_ord3(j , 1 :n) ) ; 

dist_modl ( j ,n+2) = mean(dist_modl (j , 1 :n) ) ; 



dist_mod2 ( j ,n+2) = mean(dist_mod2(j , 1 :n) ) ; 
dist_mod3(j ,n+2) = mean(dist_mod3(j , 1 :n) ) ; 

dist_orthl ( j ,n+2) = mean(dist_orthl ( j , 1 : n) ) ; 
dist_orth2(j ,n+2) = mean(dist_orth2( j , 1 : n) ) ; 
dist_orth3(j ,n+2) = mean(dist_orth3( j , 1 : n) ) ; 

dist_ordl ( j ,n+3) = var (dist_ordl ( j , 1 :n) ) ; 
dist_ord2 ( j ,n+3) = var (dist_ord2(j , 1 :n)) ; 
dist_ord3(j ,n+3) = var (dist_ord3(j , 1 :n)) ; 

dist_modl ( j ,n+3) = var (dist_modl ( j , 1 :n) ) ; 
dist_mod2 ( j ,n+3) = var (dist_mod2(j , 1 :n)) ; 
dist_mod3(j ,n+3) = var (dist_mod3(j , 1 :n)) ; 

dist_orthl ( j ,n+3) = var (dist_orthl (j , 1 :n) ) ; 
dist_orth2(j ,n+3) = var (dist_orth2 ( j , 1 :n) ) ; 
dist_orth3(j ,n+3) = var (dist_orth3(j , 1 :n) ) ; 


dist_ordl(j ,n+4) 
dist_ord2(j ,n+4) 
dist_ord3(j ,n+4) 

dist_modl(j ,n+4) 
dist_mod2(j ,n+4) 
dist_mod3(j ,n+4) 


dist_ordl ( j ,n+2) "2 
dist_ord2 ( j ,n+2) "2 
dist_ord3(j ,n+2) "2 

dist_modl ( j ,n+2) "2 
dist_mod2 ( j ,n+2) "2 
dist_mod3(j ,n+2) "2 


dist_ordl(j ,n+3) ; 
dist_ord2(j ,n+3) ; 
dist_ord3(j ,n+3) ; 

dist_modl(j ,n+3) ; 
dist_mod2(j ,n+3) ; 
dist_mod3(j ,n+3) ; 


dist_orthl(j ,n+4) 
dist_orth2(j ,n+4) 
dist_orth3(j ,n+4) 

end 


dist_orthl (j ,n+2) ~2 
dist_orth2(j ,n+2) ~2 
dist_orth3(j ,n+2) ~2 


+ dist_orthl (j ,n+3) 
+ dist_orth2(j ,n+3) 
+ dist_orth3(j ,n+3) 


distl = [dist_ordl(l :6,n+2:n+4) ; 000; 
dist_modl(l :6,n+2:n+4) ; 000; ... 
dist_orthl(l :6,n+2:n+4)] ; 
dist2 = [dist_ord2(l :6,n+2:n+4) ; 0 0 0; ... 
dist_mod2(l :6,n+2:n+4) ; 000; ... 
dist_orth2(l :6,n+2:n+4)] ; 
dist3 = [dist_ord3(l :6,n+2:n+4) ; 000; ... 
dist_mod3(l :6,n+2:n+4) ; 000; ... 
dist_orth3(l :6,n+2:n+4)] ; 


dist = [distl zeros(20,l) dist2 zeros(20,l) dist3] ; 



°/„°/o Confirmation Point Estimation 


for i = l:n 

for k = 1 : N 

for j = 1:5 

cal (i ,k) . conf ( j , : ) = ... 

mvnrnd( [xconf ( j , 1) yconf(j,l)], sigma); 
yobs (i) . conf (j ,k) = cal (i ,k) . conf (j , 2) ; 

cdist (i) . ordl (j ,k) = ( (sin(phi_ordl (j , i) )/ . 
sin(alpha_ordl ( j , i) ) ) * (yobs (i) . conf ( j ,k) 
bordl(l,i) - bordl (2 , i) *xconf ( j , 1) ) ) ; 
cdist(i) ,ord2(j ,k) = ( (sin(phi_ord2(j , i) )/ . 
sin(alpha_ord2 (j , i) ) ) * (yobs (i) . conf ( j ,k) 
bord2(l,i) - bord2(2 , i) *xconf ( j , 1) ) ) ; 
cdist(i) ,ord3(j ,k) = ( (sin(phi_ord3(j , i) )/ . 
sin(alpha_ord3(j , i) ) ) * (yobs (i) . conf ( j ,k) 
bord3(l,i) - bord3(2, i)*xconf (j , 1) ) ) ; 

cdist (i) .modi (j ,k) = ( (sin(phi_modl (j , i) )/ . 
sin(alpha_modl ( j , i) ) ) * (yobs (i) . conf ( j ,k) 
bmodl(l,i) - bmodl (2 , i) *xconf ( j , 1) ) ) ; 
cdist(i) ,mod2(j ,k) = ( (sin(phi_mod2(j , i) )/ . 
sin(alpha_mod2 (j , i) ) ) * (yobs (i) . conf ( j ,k) 
bmod2(l,i) - bmod2(2 , i) *xconf ( j , 1) ) ) ; 
cdist(i) ,mod3(j ,k) = ( (sin(phi_mod3(j , i) )/ . 
sin(alpha_mod3(j , i) ) ) * (yobs (i) . conf ( j ,k) 
bmod3(l,i) - bmod3(2 , i) *xconf ( j , 1) ) ) ; 

cdist (i) . orthl (j ,k) = ( (sin(phi_orthl ( j , i) ) / 
sin(alpha_orthl ( j , i) ) ) * (yobs (i) . conf ( j ,k) 
borthl(l,i) - borthl (2, i) *xconf ( j , 1) ) ) ; 
cdist(i) ,orth2(j ,k) = ( (sin(phi_orth2 ( j , i) ) / 
sin(alpha_orth2 ( j , i) ) ) * (yobs (i) . conf ( j ,k) 
borth2(l,i) - borth2 (2, i) *xconf ( j , 1) ) ) ; 
cdist(i) ,orth3(j ,k) = ( (sin(phi_orth3(j ,i))/ 
sin(alpha_orth3( j , i) ) ) * (yobs (i) . conf ( j ,k) 
borth3(l,i) - borth3(2, i) *xconf ( j , 1) ) ) ; 

end 


end 
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for j = 1:5 

cdist (i) . ordl ( j ,N+1) = 0; 
cdist(i) . ord2(j ,N+1) = 0; 
cdist(i) . ord3(j ,N+1) = 0; 

cdist(i) ,modl(j ,N+1) = 0; 
cdist(i) ,mod2(j ,N+1) = 0; 
cdist(i) ,mod3(j ,N+1) = 0; 

cdist (i) . orthl (j ,N+1) = 0; 
cdist (i) . orth2(j ,N+1) = 0; 
cdist (i) . orth3(j ,N+1) = 0; 

cdist (i) . ordl (j ,N+2) = raean(cdist (i) . ordl ( j , 1 : N) ) ; 
cdist(i) . ord2(j ,N+2) = raean(cdist (i) . ord2 ( j , 1 : N) ) ; 
cdist(i) . ord3(j ,N+2) = raean(cdist (i) . ord3( j , 1 : N) ) ; 

cdist (i) .modi (j ,N+2) = mean(cdist (i) .modi ( j , 1 : N) ) ; 
cdist(i) ,mod2(j ,N+2) = mean(cdist (i) ,mod2 ( j , 1 : N) ) ; 
cdist(i) ,mod3(j ,N+2) = mean(cdist (i) ,mod3( j , 1 : N) ) ; 

cdist (i) . orthl (j ,N+2) = mean(cdist (i) . orthl (j , 1 : N) ) ; 
cdist(i) . orth2(j ,N+2) = mean(cdist (i) . orth2 ( j , 1 : N) ) ; 
cdist (i) . orth3(j ,N+2) = mean(cdist (i) . orth3( j , 1 : N) ) ; 

end 

for j = 1:5 

ccdist_ordl ( j , i) = cdist (i) . ordl (j ,N+2) ; 
ccdist_ord2 ( j , i) = cdist (i) . ord2 (j ,N+2) ; 
ccdist_ord3(j , i) = cdist (i) . ord3(j ,N+2) ; 

ccdist_modl ( j , i) = cdist (i) .modi (j ,N+2) ; 
ccdist_mod2 ( j , i) = cdist(i) ,mod2(j ,N+2) ; 
ccdist_mod3(j , i) = cdist (i) ,mod3(j ,N+2) ; 

ccdist_orthl ( j , i) = cdist(i) .orthl (j ,N+2) ; 
ccdist_orth2(j , i) = cdist(i) ,orth2(j ,N+2) ; 
ccdist_orth3(j , i) = cdist(i) ,orth3(j ,N+2) ; 

end 

end 


for j = 1:5 
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ccdist_ordl (j ,n+l) = 0; 
ccdist_ord2(j ,n+l) = 0; 
ccdist_ord3(j ,n+l) = 0; 

ccdist_modl (j ,n+l) = 0; 
ccdist_mod2(j ,n+l) = 0; 
ccdist_mod3(j ,n+l) = 0; 

ccdist_orthl ( j ,n+l) = 0; 
ccdist_orth2 ( j ,n+l) = 0; 
ccdist_orth3(j ,n+l) = 0; 

ccdist_ordl (j ,n+2) = mean(ccdist_ordl ( j , 1 :n) ) ; 
ccdist_ord2(j ,n+2) = mean(ccdist_ord2(j , 1 :n) ) ; 
ccdist_ord3(j ,n+2) = mean(ccdist_ord3( j , 1 :n) ) ; 

ccdist_modl (j ,n+2) = mean(ccdist_modl ( j , 1 :n) ) ; 
ccdist_mod2(j ,n+2) = mean(ccdist_mod2(j , 1 :n) ) ; 
ccdist_mod3(j ,n+2) = mean(ccdist_mod3( j , 1 :n) ) ; 

ccdist_orthl ( j ,n+2) = mean(ccdist_orthl (j , 1 :n) ) ; 
ccdist_orth2 ( j ,n+2) = mean(ccdist_orth2 ( j , 1 :n) ) ; 
ccdist_orth3(j ,n+2) = mean(ccdist_orth3(j , 1 :n) ) ; 

ccdist_ordl (j ,n+3) = var (ccdist_ordl (j , 1 :n) ) ; 
ccdist_ord2(j ,n+3) = var (ccdist_ord2 ( j , 1 :n) ) ; 
ccdist_ord3(j ,n+3) = var(ccdist_ord3(j , 1 :n) ) ; 

ccdist_modl (j ,n+3) = var (ccdist_modl (j , 1 :n) ) ; 
ccdist_mod2(j ,n+3) = var (ccdist_mod2 ( j , 1 :n)) ; 
ccdist_mod3(j ,n+3) = var(ccdist_mod3(j , 1 :n) ) ; 

ccdist_orthl ( j ,n+3) = var (ccdist_orthl (j , 1 :n) ) ; 
ccdist_orth2 ( j ,n+3) = var (ccdist_orth2(j , 1 :n) ) ; 
ccdist_orth3(j ,n+3) = var (ccdist_orth3(j , 1 :n) ) ; 

ccdist_ordl (j ,n+4) = ccdist_ordl ( j ,n+2) "2 + ... 
ccdist_ordl ( j ,n+3) ; 

ccdist_ord2(j ,n+4) = ccdist_ord2(j ,n+2) ~2 + ... 
ccdist_ord2(j ,n+3) ; 

ccdist_ord3(j ,n+4) = ccdist_ord3(j ,n+2) ~2 + ... 
ccdist_ord3(j ,n+3) ; 

ccdist_modl (j ,n+4) = ccdist_modl ( j ,n+2) "2 + ... 
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ccdist_modl ( j ,n+3) ; 

ccdist_mod2(j ,n+4) = ccdist_mod2(j ,n+2) "2 + ... 
ccdist_mod2(j ,n+3) ; 

ccdist_mod3(j ,n+4) = ccdist_mod3(j ,n+2 ) ~2 + ... 
ccdist_mod3(j ,n+3) ; 

ccdist_orthl ( j ,n+4) = ccdist_orthl ( j ,n+2) ~2 + ... 
ccdist_orthl ( j ,n+3) ; 

ccdist_orth2 ( j ,n+4) = ccdist_orth2 ( j ,n+2) ~2 + ... 
ccdist_orth2 ( j ,n+3) ; 

ccdist_orth3(j ,n+4) = ccdist_orth3(j ,n+2) ~2 + ... 
ccdist_orth3(j ,n+3) ; 

end 

ccdistl = [ccdist_ordl (1 : 5,n+2 :n+4) ; 000; ... 
ccdist_modl (1 : 5,n+2 : n+4) ; 000; ... 
ccdist_orthl (1 :5,n+2:n+4)] ; 

ccdist2 = [ccdist_ord2(l : 5,n+2 :n+4) ; 000; ... 
ccdist_mod2(l : 5,n+2 : n+4) ; 000; ... 
ccdist_orth2 (1:5, n+2 : n+4) ] ; 

ccdist3 = [ccdist_ord3(l : 5, n+2 :n+4) ; 000; ... 
ccdist_mod3(l : 5, n+2 : n+4) ; 000; ... 
ccdist_orth3 (1:5, n+2 : n+4) ] ; 

ccdist = [ccdistl zeros(17,l) ccdist2 zeros(17,l) ccdist3] ; 

U Mean and Variance of Coefficients over n Calibrations for each 

U Estimator and Design 

bmean = [mean(bordl (1 , : ) ) mean(bmodl (1 , : ) ) mean(borthl (1 , : ) ) 0 ... 
mean(bord2 (1 , : ) ) mean(bmod2(l , : ) ) mean(borth2 (1 , : ) ) 0 ... 
mean(bord3(l , : ) ) mean(bmod3(l , : ) ) mean(borth3(l , : ) ) ; ... 
mean(bordl (2 , : ) ) mean(bmodl (2 , : ) ) mean(borthl (2 , : ) ) 0 ... 
mean(bord2 (2 , : ) ) mean(bmod2(2, : )) mean(borth2 (2 , : ) ) 0 ... 
mean(bord3(2, : )) mean(bmod3(2 , : ) ) mean(borth3 (2 ,:));] ; 

bvar = [var (bordl (1 , : ) ) var (bmodl (1 , : ) ) var (borthl (1 , : ) ) 0 ... 
var(bord2(l , : ) ) var(bmod2(l , : ) ) var (borth2(l , : ) ) 0 ... 
var(bord3(l , : ) ) var(bmod3(l , : ) ) var (borth3(l , : ) ) ; ... 
var (bordl (2 , : ) ) var (bmodl (2 , : ) ) var (borthl (2, :) ) 0 ... 
var (bord2(2 , : ) ) var(bmod2(2, : )) var (borth2(2, : ) ) 0 ... 
var (bord3(2 , : ) ) var(bmod3(2, : )) var (borth3(2, : ) ) ; ] ; 

t = toe; 

tmin = t/60; 



117 


disp( [’Elapsed time is num2str (tmin) , 


’ minutes . ’] ) ; 
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APPENDIX B 

SIMPLE QUADRATIC SIMULATION CODE 


function [dist, ccdist, bmean, bvar] = simpquadsim5(n,N,std, gamma) 

7. n: number of calibrations 

7. N: number of uses per calibration 

% std: measurement system accuracy 

7. gamma: variance ratio (delta-to-epsilon) 

warning off all; 

options = optimsetp Algorithm’ , ’ levenberg-marquardt J , . . . 

’TolFun’ ,le-12, ’TolX’ ,le-12, ’Display’ , 'off 5 ) ; 


tic ; 

7.7. Design Points 

xl = [-1; -1; 0; 0; 1; 1]; 
x2 = [-1; 0; 0; 0; 0; 1] ; 

xlm = [1-11; 1-11; 100; 100; 111; 111]; 
x2m = [1-11; 100; 100; 100; 100; 111]; 

7.7. Confirmation Points 

xconf = [-1; -0.5; 0; 0.5; 1]; 

7.7. True Model and Responses 
7. b0 = 0, bl = 1, bll = 0.2 

yl = xl + 0.2*(xl . ~2) ; 
y2 = x2 + 0 . 2* (x2 . ~2) ; 

yconf = xconf + 0 . 2* (xconf . "2) ; 

7.7. Error in x and y Sampled from Bivariate Normal Distribution 

sigma = [gamma* ( std" 2) 0; 0 (std~2)] ; 
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for i = l:n 

for j = 1:6 

cal(i) .desl(j , : ) = mvnrnd( [xl (j , 1) yl (j , 1)] , sigma) ; 
yobs_desl ( j , i) = cal(i) .desl(j ,2) ; 

cal(i) ,des2(j , : ) = mvnrnd( [x2(j , 1) y2(j , 1)] , sigma) ; 
yobs_des2(j , i) = cal(i) .des2(j ,2) ; 

end 

U Ordinary Least Squares Estimation 

bordl(:,i) = regress (yobs_desl (:, i) ,xlm) ; 
bord2(:,i) = regress (yobs_des2( i) ,x2m) ; 

U Ordinary Least Squares Distance 

for j = 1:6 

phi_ordl ( j , i) = (pi/2) - atan(bordl (2, i) ) ; 
phi_ord2 ( j , i) = (pi/2) - atan(bord2(2,i)) ; 

alpha_ordl (j , i) = (pi/2) + (1-gamma) *atan(bordl (2 , i) ) ; 
alpha_ord2(j , i) = (pi/2) + ( 1-gamma) *atan(bord2 (2, i) ) ; 

dist_ordl ( j , i) = ( (sin(phi_ordl ( j , i) ) / . . . 

sin(alpha_ordl ( j , i) ) ) * (yobs_desl ( j , i) - bordl(l,i) - ... 
bordl (2, i) *xl ( j , 1) - bordl (3, i) *xl (j , 1) *xl (j , 1) ) ) ; 
dist_ord2(j , i) = ( (sin(phi_ord2( j , i) ) / . . . 

sin(alpha_ord2(j , i) ) ) * (yobs_des2 ( j ,i) - bord2(l,i) - ... 
bord2 (2, i) *x2 ( j , 1) - bord2(3,i)*x2(j , l)*x2(j , 1))) ; 

end 

°/:/ 0 Orthogonal Least Squares Estimation 

borthl(:,i) = lsqnonlin(@(borthl) (sin( (pi/2) - ... 

atan(borthl (2) ) )/sin(pi/2) ) * (yobs_desl (1 : 6 , i) - ... 
borthl(l) - borthl (2) *xl (1 : 6 , 1) - borthl (3) * . . . 
xl (1 : 6 , 1) . ~2) , bordl (:, i) ,[],[] , options) ; 
borth2(:,i) = lsqnonlin(@(borth2) (sin( (pi/2) - ... 

atan(borth2(2) ) )/sin(pi/2) ) * (yobs_des2(l : 6 , i) - ... 
borth2 (1) - borth2 (2) *x2 (1 : 6 , 1) - borth2(3)*. . . 
x2 (1:6,1) .~2) ,bord2( :, i) ,[],[] , options) ; 



for j = 1:6 

phi_orthl ( j , i) = (pi/2) - atan(borthl (2 , i) ) ; 
phi_orth2(j , i) = (pi/2) - atan(borth2(2,i)) ; 

alpha_orthl ( j , i) = (pi/2) + (1-gamma) *atan(borthl (2 , i) ) ; 
alpha_orth2 ( j , i) = (pi/2) + (1-gamma) *atan(borth2 (2 , i) ) ; 

dist_orthl ( j , i) = ( (sin(phi_orthl ( j , i) ) / . . . 
sin(alpha_orthl (j , i) ) ) * (yobs_desl ( j , i) - borthl (l,i) - . 
borthl (2 , i) *xl ( j , 1) - borthl (3) *xl (j , 1) ~2) ) ; 
dist_orth2(j , i) = ( (sin(phi_orth2( j , i) ) / . . . 

sin(alpha_orth2(j ,i)))*(yobs_des2(j ,i) - borth2(l,i) 
borth2(2 , i) *x2( j , 1) - borth2(3) *x2(j , 1) ~2) ) ; 

end 

dist_orthl (7, i) = 0; 
dist_orth2(7, i) = 0; 

dist_orthl (8, i) = sum(dist_orthl (1 : 6, i) ) ; 
dist_orth2(8, i) = sum(dist_orth2(l : 6 , i) ) ; 

U Modified Least Squares 

bmodl (:,i) = lsqnonlin(@(bmodl) (sin( (pi/2) - ... 

atan(bmodl (2) ) ) /sin( (pi/2) + (1-gamma) *atan(bmodl (2) )))* . 
(yobs_desl (1 : 6 , i) - bmodl(l) - bmodl (2)*xl (1 : 6, 1) - ... 
bmodl (3)*xl(l:6,l) . ~2) ,bordl (:, i) ,[],[] , options) ; 
bmod2(:,i) = lsqnonlin(@(bmod2) (sin( (pi/2) - ... 

atan(bmod2(2) ) ) /sin( (pi/2) + (1-gamma) *atan(bmod2(2) )))* . 
(yobs_des2(l : 6 , i) - bmod2(l) - bmod2 (2) *x2 (1 : 6 , 1) - ... 
bmod2(3)*x2(l :6, 1) . ~2) ,bord2( :,!),[],[] , options) ; 

for j = 1:6 

phi_modl ( j , i) = (pi/2) - atan(bmodl (2, i) ) ; 
phi_mod2 ( j , i) = (pi/2) - atan(bmod2(2,i)) ; 

alpha_modl ( j , i) = (pi/2) + (1-gamma) *atan (bmodl (2 , i) ) ; 
alpha_mod2(j , i) = (pi/2) + ( 1-gamma) *atan(bmod2 (2, i) ) ; 

dist_modl ( j , i) = ( (sin(phi_modl (j , i) )/ . . . 

sin(alpha_modl ( j , i) ) ) * (yobs_desl ( j , i) - bmodl (l,i) - 
bmodl (2, i) *xl ( j , 1) - bmodl (3) *xl ( j , 1) ~2) ) ; 
dist_mod2(j , i) = ( (sin(phi_mod2(j , i) )/ . . . 



sin(alpha_mod2(j , i) ) ) * (yobs_des2 ( j ,i) - bmod2(l,i) - . 
braod2 (2, i) *x2 ( j , 1) - bmod2(3) *x2 ( j , 1) ~2) ) ; 

end 

end 

for j = 1:6 

dist_ordl ( j ,n+l) = 0; 
dist_ord2 ( j ,n+l) = 0; 

dist_modl ( j ,n+l) = 0; 
dist_mod2 ( j ,n+l) = 0; 

dist_orthl ( j ,n+l) = 0; 
dist_orth2(j ,n+l) = 0; 

dist_ordl ( j ,n+2) = mean(dist_ordl (j , 1 :n) ) ; 
dist_ord2 ( j ,n+2) = mean(dist_ord2(j , 1 :n) ) ; 

dist_modl ( j ,n+2) = mean(dist_modl (j , 1 :n) ) ; 
dist_mod2 ( j ,n+2) = mean(dist_mod2(j , 1 :n) ) ; 

dist_orthl ( j ,n+2) = mean(dist_orthl ( j , 1 : n) ) ; 
dist_orth2(j ,n+2) = mean(dist_orth2( j , 1 : n) ) ; 

dist_ordl ( j ,n+3) = var (dist_ordl ( j , 1 :n) ) ; 
dist_ord2 ( j ,n+3) = var (dist_ord2(j , 1 :n)) ; 

dist_raodl ( j ,n+3) = var (dist_modl ( j , 1 :n) ) ; 
dist_mod2 ( j ,n+3) = var (dist_mod2(j , 1 :n)) ; 

dist_orthl ( j ,n+3) = var (dist_orthl (j , 1 :n) ) ; 
dist_orth2(j ,n+3) = var (dist_orth2 ( j , 1 :n) ) ; 

dist_ordl ( j ,n+4) = dist_ordl ( j ,n+2) ~2 + dist_ordl (j ,n+3) ; 
dist_ord2 ( j ,n+4) = dist_ord2 ( j ,n+2) ~2 + dist_ord2(j ,n+3) ; 

dist_raodl ( j ,n+4) = dist_modl ( j ,n+2) ~2 + dist_modl (j ,n+3) ; 
dist_raod2 ( j ,n+4) = dist_mod2 ( j ,n+2) ~2 + dist_mod2(j ,n+3) ; 

dist_orthl ( j ,n+4) = dist_orthl (j ,n+2) "2 + dist_orthl ( j ,n+3) ; 
dist_orth2(j ,n+4) = dist_orth2(j ,n+2) ~2 + dist_orth2(j ,n+3) ; 


end 



distl = [dist_ordl(l :6,n+2:n+4) ; 000; ... 
dist_modl(l :6,n+2:n+4) ; 000; ... 
dist_orthl(l :6,n+2:n+4)] ; 
dist2 = [dist_ord2(l :6,n+2:n+4) ; 000; ... 
dist_mod2(l :6,n+2:n+4) ; 000; ... 
dist_orth2(l :6,n+2:n+4)] ; 

dist = [distl zeros(20,l) dist2] ; 

for i = l:n 

for k = 1 : N 

for j = 1:5 

cal (i ,k) . conf ( j , : ) = ... 

mvnrnd( [xconf ( j , 1) yconf(j,l)], sigma); 
yobs (i) . conf (j ,k) = cal (i ,k) . conf (j , 2) ; 

cdist (i) . ordl (j ,k) = ( (sin(phi_ordl (j , i) )/ . . 
sin(alpha_ordl ( j , i) ) ) * (yobs (i) . conf ( j ,k) 
bordl (l,i) - bordl (2 , i) *xconf ( j , 1) - ... 
bordl (3, i) *xconf (j , 1) ~2) ) ; 
cdist(i) ,ord2(j ,k) = ( (sin(phi_ord2(j , i) )/ . . 
sin(alpha_ord2 (j , i) ) ) * (yobs (i) . conf ( j ,k) 
bord2(l,i) - bord2(2 , i) *xconf ( j , 1) - ... 
bord2(3, i) *xconf (j , 1) ~2) ) ; 

cdist (i) .modi (j ,k) = ( (sin(phi_modl (j , i) )/ . . 
sin(alpha_modl (j , i) ) ) * (yobs (i) . conf ( j ,k) 
bmodl (l,i) - bmodl (2 , i) *xconf ( j , 1) - ... 
bmodl (3, i) *xconf (j , 1) ~2) ) ; 
cdist(i) ,mod2(j ,k) = ( (sin(phi_mod2(j , i) )/ . . 
sin(alpha_mod2 ( j , i) ) ) * (yobs (i) . conf ( j ,k) 
bmod2(l,i) - bmod2(2 , i) *xconf ( j , 1) - ... 
bmod2(3, i) *xconf (j , 1) ~2) ) ; 

cdist (i) . orthl (j ,k) = ( (sin(phi_orthl ( j , i) ) / 
sin(alpha_orthl ( j , i) ) ) * (yobs (i) . conf ( j ,k) 
borthl(l,i) - borthl (2, i) *xconf ( j , 1) - .. 
borthl (3 , i) *xconf ( j , 1) ~2) ) ; 
cdist(i) ,orth2(j ,k) = ( (sin(phi_orth2 ( j , i) ) / 
sin(alpha_orth2 ( j , i) ) ) * (yobs (i) . conf ( j ,k) 
borth2(l,i) - borth2 (2, i) *xconf ( j , 1) - .. 
borth2 (3 , i) *xconf ( j , 1) ~2) ) ; 
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end 

end 

for j = 1:5 

cdist(i) .ordl(j ,N+1) = 0; 
cdist(i) . ord2(j ,N+1) = 0; 

cdist (i) .modi (j ,N+1) = 0; 
cdist(i) ,mod2(j ,N+1) = 0; 

cdist (i) . orthl (j ,N+1) = 0; 
cdist (i) . orth2(j ,N+1) = 0; 

cdist (i) . ordl (j ,N+2) = mean(cdist (i) . ordl ( j , 1 : N) ) ; 
cdist(i) . ord2(j ,N+2) = mean(cdist (i) . ord2 ( j , 1 : N) ) ; 

cdist (i) .modi (j ,N+2) = mean(cdist (i) .modi (j , 1 :N) ) ; 
cdist(i) ,mod2(j ,N+2) = mean(cdist (i) ,mod2 ( j , 1 : N) ) ; 

cdist(i) . orthl(j ,N+2) = mean(cdist (i) . orthl (j , 1 : N) ) ; 
cdist(i) . orth2(j ,N+2) = mean(cdist (i) . orth2 ( j , 1 : N) ) ; 

cdist (i) . ordl (j ,N+3) = var (cdist (i) . ordl (j , 1 : N) ) ; 
cdist(i) . ord2(j ,N+3) = var (cdist (i) . ord2(j , 1 : N) ) ; 

cdist(i) ,modl(j ,N+3) = var (cdist (i) .modi (j , 1 : N) ) ; 
cdist(i) ,mod2(j ,N+3) = var (cdist (i) ,mod2(j , 1 : N) ) ; 

cdist (i) . orthl (j ,N+3) = var (cdist (i) . orthl (j , 1 : N) ) ; 
cdist (i) . orth2(j ,N+3) = var (cdist (i) . orth2(j , 1 : N) ) ; 

end 

for j = 1:5 

ccdist_ordl ( j , i) = cdist (i) . ordl (j ,N+2) ; 
ccdist_ord2 ( j , i) = cdist (i) . ord2 (j ,N+2) ; 

ccdist_modl ( j , i) = cdist (i) .modi (j ,N+2) ; 
ccdist_mod2 ( j , i) = cdist(i) ,mod2(j ,N+2) ; 

ccdist_orthl ( j , i) = cdist(i) .orthl (j ,N+2) ; 
ccdist_orth2(j , i) = cdist (i) ,orth2(j ,N+2) ; 


end 
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end 

for j = 1:5 

ccdist_ordl (j ,n+l) = 0; 
ccdist_ord2(j ,n+l) = 0; 

ccdist_modl (j ,n+l) = 0; 
ccdist_mod2(j ,n+l) = 0; 

ccdist_orthl ( j ,n+l) = 0; 
ccdist_orth2 ( j ,n+l) = 0; 

ccdist_ordl (j ,n+2) = mean(ccdist_ordl ( j , 1 :n) ) ; 
ccdist_ord2(j ,n+2) = mean(ccdist_ord2(j , 1 :n) ) ; 

ccdist_modl (j ,n+2) = mean(ccdist_modl ( j , 1 :n) ) ; 
ccdist_mod2(j ,n+2) = mean(ccdist_mod2(j , 1 :n) ) ; 

ccdist_orthl ( j ,n+2) = mean(ccdist_orthl (j , 1 :n) ) ; 
ccdist_orth2 ( j ,n+2) = mean(ccdist_orth2 ( j , 1 :n)) ; 

ccdist_ordl (j ,n+3) = var (ccdist_ordl (j , 1 :n) ) ; 
ccdist_ord2(j ,n+3) = var (ccdist_ord2 ( j , 1 :n) ) ; 

ccdist_modl (j ,n+3) = var (ccdist_modl (j , 1 :n) ) ; 
ccdist_mod2(j ,n+3) = var (ccdist_mod2 ( j , 1 :n)) ; 

ccdist_orthl ( j ,n+3) = var (ccdist_orthl (j , 1 :n) ) ; 
ccdist_orth2 ( j ,n+3) = var (ccdist_orth2(j , 1 :n) ) ; 

ccdist_ordl (j ,n+4) = ccdist_ordl ( j ,n+2) "2 + ccdist_ordl ( j ,n+3) ; 
ccdist_ord2(j ,n+4) = ccdist_ord2(j ,n+2) ~2 + ccdist_ord2 ( j ,n+3) ; 

ccdist_modl (j ,n+4) = ccdist_modl ( j ,n+2 ) ~2 + ccdist_modl ( j ,n+3) ; 
ccdist_mod2(j ,n+4) = ccdist_mod2(j ,n+2) ~2 + ccdist_mod2 ( j ,n+3) ; 

ccdist_orthl ( j ,n+4) = ccdist_orthl ( j ,n+2) ~2 + ccdist_orthl ( j ,n+3) ; 
ccdist_orth2 ( j ,n+4) = ccdist_orth2 ( j ,n+2) ~2 + ccdist_orth2 ( j ,n+3) ; 

end 

ccdistl = [ccdist_ordl (1 : 5,n+2 : n+4) ; 000; ... 
ccdist_modl (1 : 5,n+2 : n+4) ; 000; ... 
ccdist_orthl (1 :5,n+2:n+4)] ; 



ccdist2 = [ccdist_ord2(l : 5,n+2 :n+4) ; 000; ... 
ccdist_mod2(l : 5,n+2 : n+4) ; 000; ... 
ccdist_orth2(l :5,n+2:n+4)] ; 

ccdist = [ccdistl zeros(17,l) ccdist2] ; 

°/„°/o Mean and Variance of Coefficients over n Calibrations for 
°/„°/o each Estimator and Design 

bmean = [mean(bordl (1 , : ) ) mean(bmodl (1 , : ) ) mean(borthl (1 , : ) ) 0 
mean(bord2 (1 , : ) ) mean(bmod2(l , : ) ) mean(borth2 (1 , : ) ) ; ... 
mean(bordl (2 , : ) ) mean(bmodl (2 , : ) ) mean(borthl (2 , : ) ) 0 ... 
mean(bord2 (2 , : ) ) mean(bmod2(2, : )) mean(borth2 (2 ,:));] ; 
bvar = [var (bordl (1 , : ) ) var (bmodl (1 , : ) ) var (borthl (1 ,:) ) 0 ... 
var(bord2(l , : ) ) var(bmod2(l , : ) ) var (borth2(l , : ) ) ; ... 
var (bordl (2 , : ) ) var (bmodl (2 , : ) ) var (borthl (2, :) ) 0 ... 
var (bord2(2 , : ) ) var(bmod2(2, : )) var (borth2(2, : ) ) ; ] ; 

t = toe; 
tmin = t/60; 

disp( ['Elapsed time is num2str (tmin) , ' minutes.']); 



126 


APPENDIX C 


6-COMPONENT FORCE-BALANCE SIMULATION 

CODE 


function [dist, ccdist, bmean, bvar] = mult sim5(n,N,std, gamma) 

°/„ n: number of calibrations 

°/„ N: number of uses per calibration 

°/ 0 std: measurement system accuracy 

°/ 0 gamma: variance ratio (delta-to-epsilon) 

warning off all; 

options = optimset(’ Algorithm’ , ’ levenberg-marquardt ’ , . . . 

’ TolFun’ , le-12, ’TolX’ , le-12, ’Display’ , 'off 5 ) ; 


tic ; 


U Design Points - Central Composite Design 


xl = [-0.639 -0.639 -0.639 -0.639 -0.639 -0.639; ... 
0.000 0.000 0.000 0.000 0.000 0.000; ... 

-0.639 -0.639 0.639 0.639 0.639 0.639; ... 

0.000 -1.00 0.000 0.000 0.000 0.000; ... 

0.639 0.639 0.639 -0.639 -0.639 0.639; ... 

-0.639 0.639 -0.639 0.639 0.639 0.639; ... 

0.639 0.639 0.639 0.639 0.639 0.639; ... 

0.639 -0.639 -0.639 0.639 -0.639 -0.639; ... 

0.000 0.000 0.000 0.000 0.000 1.00; ... 

-0.639 -0.639 0.639 0.639 -0.639 -0.639; ... 

0.639 -0.639 0.639 -0.639 -0.639 -0.639; ... 

0.000 0.000 0.000 0.000 0.000 0 . 000 ; ... 

1.00 0.000 0.000 0.000 0.000 0.000; ... 

0.639 0.639 -0.639 -0.639 0.639 0.639; ... 

0.639 -0.639 -0.639 -0.639 -0.639 0.639; ... 

0.639 -0.639 0.639 0.639 -0.639 0.639; ... 

0.000 0.000 0.000 0.000 0.000 -1.00; ... 

-0.639 0.639 0.639 0.639 0.639 -0.639; ... 
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- 1.00 0.000 0.000 0.000 0.000 0 . 000 ; ... 
0.000 0.000 0.000 0.000 0.000 0 . 000 ; ... 
- 0.639 - 0.639 - 0.639 0.639 - 0.639 0 . 639 ; 
0.639 0.639 - 0.639 - 0.639 - 0.639 - 0 . 639 ; 
0.639 0.639 0.639 - 0.639 0.639 - 0 . 639 ; . 
- 0.639 0.639 0.639 0.639 - 0.639 0 . 639 ; . 
- 0.639 0.639 - 0.639 0.639 - 0.639 - 0 . 639 ; 
0.000 0.000 0.000 - 1.00 0.000 0 . 000 ; ... 

0.000 0.000 0.000 0.000 0.000 0 . 000 ; ... 

- 0.639 0.639 0.639 - 0.639 - 0.639 - 0 . 639 ; 
0.639 - 0.639 - 0.639 - 0.639 0.639 - 0 . 639 ; 
0.000 1.00 0.000 0.000 0.000 0 . 000 ; ... 
- 0.639 - 0.639 - 0.639 - 0.639 0.639 0 . 639 ; 
0.639 0.639 - 0.639 0.639 0.639 - 0 . 639 ; . 
0.639 - 0.639 0.639 0.639 0.639 - 0 . 639 ; . 
0.639 0.639 - 0.639 0.639 - 0.639 0 . 639 ; . 
0.639 - 0.639 0.639 - 0.639 0.639 0 . 639 ; . 
0.000 0.000 0.000 0.000 0.000 0 . 000 ; ... 
- 0.639 - 0.639 0.639 - 0.639 0.639 - 0 . 639 ; 
0.000 0.000 0.000 0.000 0.000 0 . 000 ; ... 

0.000 0.000 0.000 0.000 1.00 0 . 000 ; ... 

0.000 0.000 0.000 1.00 0.000 0 . 000 ; ... 

0.639 - 0.639 - 0.639 0.639 0.639 0 . 639 ; . 
0.000 0.000 0.000 0.000 0.000 0 . 000 ; ... 

0.000 0.000 0.000 0.000 - 1.00 0 . 000 ; ... 

- 0.639 0.639 0.639 - 0.639 0.639 0 . 639 ; . 
0.639 0.639 0.639 0.639 - 0.639 - 0 . 639 ; . 
0.000 0.000 0.000 0.000 0.000 0 . 000 ; ... 
- 0.639 0.639 - 0.639 - 0.639 0.639 - 0 . 639 ; 
0.000 0.000 1.00 0.000 0.000 0 . 000 ; ... 
0.000 0.000 0.000 0.000 0.000 0 . 000 ; ... 
- 0.639 - 0.639 0.639 - 0.639 - 0.639 0 . 639 ; 
- 0.639 - 0.639 - 0.639 0.639 0.639 - 0 . 639 ; 
0.000 0.000 0.000 0.000 0.000 0 . 000 ; ... 
- 0.639 0.639 - 0.639 - 0.639 - 0.639 0 . 639 ; 
0.000 0.000 - 1.00 0.000 0.000 0 . 000 ]; 


xlm = x 2 fx ( xl , ’ quadratic ’ ) ; 
°/„°/o Confirmation Points 


xconf = [- 0.138 - 0.600 - 0.414 0.264 - 0.138 - 0 . 600 ; ... 
- 0.600 - 0.600 - 0.600 0.600 0.600 0 . 366 ; ... 

- 0.600 - 0.600 0.372 - 0.600 - 0.600 - 0 . 600 ; ... 



-0.126 0.600 -0.0960 -0.600 -0.228 -0.270; ... 
-0.126 0.600 -0.0960 -0.600 -0.228 -0.270; ... 
-0.600 -0.600 0.600 0.600 0.600-0.600; ... 
-0.600 0.600 0.600 -0.600 -0.600 0.600; ... 
-0.528 0.600 -0.348 0.288 -0.00600 0.600; ... 
0.600 -0.600 -0.600 -0.600 -0.258 0.600; ... 
-0.600 0.600 -0.600 0.600 -0.600 -0.600; ... 
0.600 0.228 -0.552 0.600 0.600 0.600; ... 
-0.540 -0.408 0.0720 -0.474 0.462 0.600; ... 
-0.450 -0.120 -0.600 -0.288 -0.600 0.168; ... 
-0.600 -0.124 0.150 0.600 -0.264 -0.0480; ... 
-0.186 0.174 0.600 0.000 -0.306 -0.600; ... 
-0.600 -0.124 0.150 0.600 -0.264 -0.0480; ... 
0.504 -0.558 -0.0942 -0.0780 0.600 -0.0840; .. 
0.600 0.600 0.600 -0.600 0.600 -0.600; ... 
0.600 -0.120 -0.600 -0.600 -0.600 -0.600; ... 
0.600 0.528 -0.600 0.234 0.0809 -0.216; ... 
0.600 0.600 0.162 0.600 -0.227 -0.600; ... 
0.0308 0.252 -0.0715 0.600 0.600-0.510; ... 
0.180 -0.270 0.600 -0.600 0.102 0.0360; ... 
-0.126 -0.600 0.600 0.239 -0.600 0.600; ... 
0.600 0.216 0.178 -0.228 -0.600 0.456; ... 
0.180 -0.270 0.600 -0.600 0.102 0.0360; ... 
0.288 0.540 -0.600 -0.600 0.600 0.600; ... 
0.600 -0.600 -0.324 0.0420 -0.600 -0.0180; ... 
0.504 -0.558 -0.0942 -0.0780 0.600 -0.0840; .. 
0.600 -0.600 0.600 0.600 -0.600 -0.534; ... 
0.366 -0.180 -0.600 0.600 -0.600 0.600; ... 
-0.186 0.174 0.600 0.000 -0.306 -0.600; ... 
-0.576 0.600 0.582 -0.0660 0.600 0.0720; ... 
0.240 0.600 0.600 0.600 -0.600 0.240; ... 
-0.600 0.540 -0.600 -0.466 0.114 0.250; ... 
-0.600 0.0420 -0.600 -0.600 0.600-0.600; ... 
0.600 0.600 0.126 -0.0900 0.336 0.341; ... 
0.600 -0.0360 0.600 0.600 0.402 0.600]; 


xconfm = x2fx(xconf , ’quadratic’ ) ; 

°/„°/o True Model and Responses 

7. Beta's Based on NTF-113C Calibration - Axial Force Component 


beta = [0 0.1334 1 -0.0095 0.0257 -0.0029 0.0006 0.0321 ... 
0.0934 0.0014 -0.0027 0.0023 -0.0025 0 -0.0003 0.0009 0 ... 
-0.0013 -0.0001 -0.0683 -0.0105 0.0306 -0.1015 0.0007 0.0387 
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0.0624 -0.0025 -0.0299] ’ ; 
yl = xlm*beta; 
yconf = xconfm*beta; 

U Error in x and y Sampled from Multi-variate 

°/„°/o Normal Distribution 

sigma = [gamma* (std~2) 000000; 0 gamma* (std~2) 00000; 

0 0 gamma* (std~2) 0000; 000 gamma* (std~2) 000; ... 

0000 gamma* (std~2) 00; 00000 gamma* (std~2) 0; ... 

000000 std~2] ; 

for i = l:n 

for j = 1 : length(xl ( : , 1) ) 

cal(i) ,desl(j , : ) = mvnrnd( [xl (j , 1) xl(j,2) xl(j,3) ... 

xl ( j , 4) xl ( j , 5) xl ( j , 6) yl ( j ,1)] , sigma) ; 
yobs_desl ( j , i) = cal(i) ,desl(j ,7) ; 

end 

°/:/ 0 Ordinary Least Squares Estimation 

bordl(:,i) = regress (yobs_desl (:, i) ,xlm) ; 

°/:/ 0 Ordinary Least Squares Distance 

for j = 1 : length(xl ( : , 1) ) 

phi_ordl ( j , i) = (pi/2) - atan(sqrt( ... 

(bordl (2 , i) + bordl (8, i) *xl ( j , 2) + ... 
bordl (9 , i) *xl ( j ,3) + bordl (10 , i) *xl ( j ,4) + ... 
bordl (11 ,i)*xl(j,5) + bordl(12,i)*xl(j,6) + ... 
2*bordl (23,i)*xl(j,l))~2 + (bordl (3, i) + ... 
bordl(8,i)*xl(j , 1) + bordl (13, i) *xl ( j ,3) + ... 
bordl(14,i)*xl(j ,4) + bordl(15,i)*xl(j ,5) + ... 
bordl (16 ,i)*xl(j,6) + 2*bordl (24, i) *xl ( j , 2) ) ~2 + . 
(bordl (4, i) + bordl (9 , i) *xl ( j , 1) + ... 
bordl (13 ,i)*xl(j,2) + bordl(17,i)*xl(j ,4) + ... 
bordl (18 , i) *xl (j , 5) + bordl(19,i)*xl(j ,6) + ... 
2*bordl (25 , i) *xl ( j ,3) ) ~2 + (bordl(5,i) + ... 
bordl (10 , i) *xl (j , 1) + bordl (14, i) *xl (j , 2) + ... 
bordl(17,i)*xl(j ,3) + bordl(20,i)*xl(j ,5) + ... 
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bordl (21,i)*xl(j,6) + 2*bordl (26 , i) *xl ( j ,4) ) ~2 + ... 
(bordl(6,i) + bordl (11 , i) *xl (j , 1) + ... 
bordl(15,i)*xl(j ,2) + bordl (18 , i) *xl ( j , 3) + ... 
bordl(20,i)*xl(j ,4) + bordl (22 , i) *xl ( j , 6) + ... 

2*bordl (27,i)*xl(j,5))~2 + (bordl (7, i) + ... 
ordl (12 , i) *xl ( j , 1) + bordl (16 , i) *xl ( j ,2) + ... 
bordl(19,i)*xl(j ,3) + bordl (21 , i) *xl (j ,4) + ... 
bordl(22,i)*xl(j ,5) + 2*bordl (28 , i) *xl ( j , 6) ) ~2) ) ; 

alpha_ordl ( j , i) = (pi/2) + (1-gamma) *atan(sqrt ( ... 

(bordl (2, i) + bordl (8, i) *xl ( j , 2) + ... 

bordl (9 , i) *xl ( j ,3) + bordl (10 , i) *xl ( j ,4) + ... 

bordl (11 , i)*xl(j,5) + bordl (12 , i) *xl ( j , 6) + ... 

2*bordl (23,i)*xl(j,l))~2 + (bordl (3, i) + ... 
bordl (8, i) *xl ( j , 1) + bordl (13, i) *xl ( j ,3) + ... 
bordl (14, i) *xl (j ,4) + bordl(15,i)*xl(j ,5) + ... 
bordl (16, i)*xl(j ,6) + 2*bordl (24, i) *xl ( j , 2) ) ~2 + ... 
(bordl (4, i) + bordl (9 , i) *xl (j , 1) + ... 
bordl(13,i)*xl(j ,2) + bordl(17,i)*xl(j ,4) + ... 
bordl(18,i)*xl(j ,5) + bordl (19 , i) *xl ( j , 6) + ... 

2*bordl (25 , i) *xl ( j ,3) ) "2 + (bordl(5,i) + ... 
bordl (10, i)*xl(j , 1) + bordl(14,i)*xl(j ,2) + ... 
bordl(17,i)*xl(j ,3) + bordl (20 , i) *xl ( j , 5) + ... 
bordl (21 , i) *xl (j , 6) + 2*bordl (26 , i) *xl ( j ,4) ) "2 + ... 
(bordl(6,i) + bordl (11 , i) *xl (j , 1) + ... 
bordl(15,i)*xl(j ,2) + bordl (18 , i) *xl ( j , 3) + ... 
bordl(20,i)*xl(j ,4) + bordl (22 , i) *xl ( j , 6) + ... 

2*bordl (27,i)*xl(j,5))"2 + (bordl (7, i) + ... 
bordl(12,i)*xl(j , 1) + bordl (16 , i) *xl ( j , 2) + ... 
bordl (19, i)*xl(j ,3) + bordl (21, i)*xl(j ,4) + ... 
bordl (22, i)*xl(j ,5) + 2*bordl (28 , i) *xl ( j , 6) ) ~2) ) ; 

ratio_ordl ( j , i) = sin(phi_ordl ( j , i) ) /sin(alpha_ordl ( j , i) ) ; 


end 

res_ordl (1 : j , i) = (yobs_desl ( : , i) - xlm*bordl ( : , i) ) ; 
for j = 1 : length(xl ( : , 1) ) 

dist_ordl ( j , i) = ratio_ordl ( j , i) *res_ordl ( j , i) ; 


end 
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U Orthogonal Least Squares Estimation 

borthl(:,i) = lsqnonlin(@(borthl) (sin( (pi/2) - atan(sqrt( ... 
(borthl (2) + borthl (8)*xl(l:j,2) + borthl (9) *xl (1 : j , 3) + ... 
borthl(10)*xl(l: j ,4) + borthl (11) *xl (1 : j , 5) + ... 
borthl (12) *xl (1 : j ,6) + 2*borthl (23) *xl (1 : j , 1) ) . “2 + ... 
(borthl(3) + borthl (8) *xl (1 : j , 1) + borthl (13) *xl (1 : j ,3) + ... 
borthl (14) *xl (1 : j ,4) + borthl (15) *xl (1 : j , 5) + ... 
borthl (16) *xl (1 : j , 6) + 2*borthl (24) *xl (1 : j , 2) ) . ~2 + ... 

(borthl (4) + borthl (9)*xl(l:j,l) + borthl (13) *xl (1 : j ,2) + ... 
borthl (17) *xl (1 : j ,4) + borthl (18) *xl (1 : j , 5) + ... 
borthl (19) *xl (1 : j ,6) + 2*borthl(25)*xl(l : j ,3)) . “2 + ... 

(borthl (5) + borthl(10)*xl(l:j,l) + borthl (14) *xl (1 : j , 2) + ... 
borthl (17) *xl (1 : j ,3) + borthl (20) *xl (1 : j , 5) + ... 
borthl (21)*xl(l:j,6) + 2*borthl (26) *xl (1 : j ,4)) . “2 + ... 
(borthl(6) + borthl (11) *xl (1 : j , 1) + borthl (15) *xl (1 : j , 2) + ... 
borthl (18) *xl (1 : j ,3) + borthl (20) *xl (1 : j ,4) + ... 
borthl (22)*xl(l:j,6) + 2*borthl (27)*xl (1 : j ,5) ) . “2 + ... 

(borthl (7) + borthl (12) *xl (1 : j , 1) + borthl (16) *xl (1 : j , 2) + ... 
borthl (19) *xl (1 : j ,3) + borthl (21) *xl (1 : j ,4) + ... 
borthl (22)*xl(l:j,5) + 2*borthl (28) *xl (1 : j , 6) ) . “2) ) ) . / . . . 
sin(pi/2) ) . * (yobs_desl (1 : j , i) - xlm*borthl (1 : 28) ) , . . . 
bordl (:, i) ,[],[] , options) ; 

for j = 1 : length(xl ( : , 1) ) 

phi_orthl ( j , i) = (pi/2) - atan(sqrt( ... 

(borthl(2,i) + borthl (8 , i) *xl ( j , 2) + ... 
borthl (9 , i) *xl (j ,3) + borthl (10 , i) *xl (j ,4) + ... 
borthl (11 , i) *xl (j , 5) + borthl (12 , i) *xl ( j , 6) + ... 

2*borthl (23,i)*xl(j,l))~2 + (borthl (3, i) + ... 
borthl(8,i)*xl(j , 1) + borthl (13 , i) *xl ( j , 3) + ... 
borthl (14, i)*xl (j ,4) + borthl(15,i)*xl(j ,5) + ... 
borthl (16, i)*xl(j ,6) + 2*borthl (24, i) *xl ( j , 2) ) "2 + ... 
(borthl (4, i) + borthl (9 , i) *xl (j , 1) + ... 
borthl (13, i)*xl(j ,2) + borthl(17,i)*xl(j ,4) + ... 
borthl(18,i)*xl(j ,5) + borthl (19 , i) *xl ( j , 6) + ... 

2*borthl (25 , i) *xl ( j ,3) ) ~2 + (borthl(5,i) + ... 
borthl (10, i) *xl (j , 1) + borthl (14, i) *xl ( j ,2) + ... 
borthl(17,i)*xl(j ,3) + borthl (20 , i) *xl (j ,5) + ... 
borthl (21,i)*xl(j,6) + 2*borthl (26, i) *xl ( j ,4) ) "2 + ... 
(borthl(6,i) + borthl (11 , i) *xl (j , 1) + ... 
borthl(15,i)*xl(j ,2) + borthl(18,i)*xl(j ,3) + ... 
borthl (20, i)*xl(j ,4) + borthl (22, i)*xl(j ,6) + ... 

2*borthl (27,i)*xl(j,5))~2 + (borthl (7, i) + ... 



borthl(12,i)*xl(j , 1) + borthl (16 , i) *xl ( j ,2) + ... 
borthl(19,i)*xl(j ,3) + borthl (21 , i) *xl ( j ,4) + ... 
borthl (22, i)*xl(j ,5) + 2*borthl (28, i) *xl ( j , 6) ) ~2) ) 

alpha_orthl ( j , i) = (pi/2) + (1-gamma) *atan(sqrt ( ... 
(borthl(2,i) + borthl(8,i)*xl(j ,2) + ... 
borthl (9 , i) *xl (j ,3) + borthl (10 , i) *xl (j ,4) + ... 
borthl (11 , i) *xl (j , 5) + borthl (12 , i) *xl (j ,6) + ... 
2*borthl (23, i) *xl (j , 1) ) ~2 + (borthl(3,i) + ... 
borthl (8, i)*xl (j , 1) + borthl (13 , i) *xl ( j , 3) + ... 
borthl (14, i) *xl (j ,4) + borthl(15,i)*xl(j ,5) + ... 
borthl (16, i)*xl(j ,6) + 2*borthl (24, i) *xl ( j , 2) ) "2 + 
(borthl(4,i) + borthl(9,i)*xl(j , 1) + ... 
borthl (13, i)*xl(j ,2) + borthl(17,i)*xl(j ,4) + ... 
borthl(18,i)*xl(j ,5) + borthl (19, i) *xl(j ,6) + ... 
2*borthl (25 , i) *xl ( j ,3) ) ~2 + (borthl(5,i) + ... 
borthl (10, i)*xl(j , 1) + borthl (14, i) *xl ( j ,2) + ... 
borthl(17,i)*xl(j ,3) + borthl (20 , i) *xl (j ,5) + ... 
borthl (21,i)*xl(j,6) + 2*borthl (26, i) *xl ( j ,4) ) "2 + 
(borthl(6,i) + borthl (11 , i) *xl (j , 1) + ... 
borthl(15,i)*xl(j ,2) + borthl(18,i)*xl(j ,3) + ... 
borthl (20, i)*xl(j ,4) + borthl (22, i)*xl(j ,6) + ... 
2*borthl (27,i)*xl(j,5))~2 + (borthl (7, i) + ... 
borthl(12,i)*xl(j , 1) + borthl (16 , i) *xl (j ,2) + ... 
borthl (19, i)*xl(j ,3) + borthl (21 , i) *xl ( j ,4) + ... 
borthl (22, i)*xl(j ,5) + 2*borthl (28, i) *xl ( j , 6) ) “2) ) 

ratio_orthl ( j , i) = sin(phi_orthl ( j , i) ) / . . . 
sin(alpha_orthl (j , i) ) ; 


end 

res_orthl (1 : j , i) = (yobs_desl ( : , i) - xlm*borthl ( : , i) ) ; 
for j = 1 : length(xl ( : , 1) ) 

dist_orthl (j , i) = ratio_orthl (j , i) *res_orthl ( j , i) ; 

end 

U Modified Least Squares 

bmodl(:,i) = lsqnonlin(@(bmodl) (sin( (pi/2) - atan(sqrt( . 
(bmodl(2) + bmodl (8) *xl (1 : j , 2) + bmodl (9) *xl (1 : j , 3) + 
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bmodl (10)*xl (1 : j ,4) + bmodl (11) *xl (1 : j , 5) + ... 
bmodl (12) *xl (1 : j , 6) + 2*bmodl (23) *xl (1 : j , 1) ) . ~2 + ... 

(bmodl (3) + bmodl (8) *xl (1 : j , 1) + bmodl (13) *xl (1 : j ,3) + ... 
bmodl (14) *xl (1 : j ,4) + bmodl (15) *xl (1 : j , 5) + ... 
bmodl (16) *xl (1 : j , 6) + 2*bmodl (24) *xl (1 : j , 2) ) . ~2 + ... 

(bmodl (4) + bmodl (9) *xl (1 : j , 1) + bmodl (13) *xl (1 : j , 2) + ... 
bmodl (17) *xl (1 : j ,4) + bmodl (18)*xl (1 : j ,5) + ... 
bmodl (19) *xl (1 : j , 6) + 2*bmodl (25) *xl (1 : j , 3) ) . “2 + ... 

(bmodl (5) + bmodl (10)*xl (1 : j , 1) + bmodl (14)*xl (1 : j ,2) + ... 
bmodl (17)*xl (1 : j ,3) + bmodl (20) *xl (1 : j , 5) + ... 
bmodl (21)*xl (1 : j ,6) + 2*bmodl (26) *xl (1 : j ,4)) . ~2 + ... 

(bmodl (6) + bmodl (ll)*xl (1 : j , 1) + bmodl (15)*xl (1 : j ,2) + ... 
bmodl (18)*xl (1 : j ,3) + bmodl (20) *xl (1 : j ,4) + ... 
bmodl (22) *xl(l: j ,6) + 2*bmodl (27) *xl (1 : j , 5) ) . “2 + ... 

(bmodl (7) + bmodl (12) *xl (1 : j , 1) + bmodl (16)*xl (1 : j ,2) + ... 
bmodl (19) *xl (1 : j ,3) + bmodl (21) *xl (1 : j ,4) + ... 
bmodl (22) *xl (1 : j , 5) + 2*bmodl (28) *xl (1 : j , 6) ) . ~2 ... 

) ) ) . /sin( (pi/2) + (1-gamma) *atan(sqrt ( ... 

(bmodl (2) + bmodl (8) *xl (1 : j ,2) + bmodl (9) *xl (1 : j , 3) + ... 
bmodl (10)*xl (1 : j ,4) + bmodl (11) *xl (1 : j , 5) + ... 
bmodl (12) *xl (1 : j , 6) + 2*bmodl (23) *xl (1 : j , 1) ) . “2 + ... 

(bmodl (3) + bmodl (8) *xl (1 : j , 1) + bmodl (13) *xl (1 : j ,3) + ... 
bmodl (14)*xl (1 : j ,4) + bmodl (15) *xl (1 : j , 5) + ... 
bmodl (16)*xl (1 : j , 6) + 2*bmodl (24) *xl (1 : j , 2) ) . ~2 + ... 

(bmodl (4) + bmodl (9) *xl (1 : j , 1) + bmodl (13) *xl (1 : j , 2) + ... 
bmodl (17) *xl (1 : j ,4) + bmodl (18) *xl (1 : j , 5) + ... 
bmodl (19) *xl (1 : j , 6) + 2*bmodl (25) *xl (1 : j , 3) ) . “2 + ... 

(bmodl (5) + bmodl (10)*xl (1 : j , 1) + bmodl (14)*xl (1 : j ,2) + ... 
bmodl (17)*xl (1 : j ,3) + bmodl (20) *xl (1 : j , 5) + ... 
bmodl (21)*xl (1 : j , 6) + 2*bmodl (26) *xl (1 : j ,4)) . ~2 + ... 

(bmodl (6) + bmodl (11) *xl (1 : j , 1) + bmodl (15)*xl (1 : j ,2) + ... 
bmodl (18)*xl (1 : j ,3) + bmodl (20) *xl (1 : j ,4) + ... 
bmodl (22)*xl(l:j,6) + 2*bmodl (27) *xl (1 : j , 5) ) . ~2 + ... 

(bmodl (7) + bmodl (12)*xl (1 : j , 1) + bmodl (16) *xl (1 : j ,2) + ... 
bmodl (19) *xl (1 : j ,3) + bmodl (21) *xl (1 : j ,4) + ... 
bmodl (22) *xl (1 : j , 5) + 2*bmodl (28) *xl (1 : j , 6) ) . “2 ... 

)))).* (yobs_desl (1 : j , i) - xlm*bmodl (1 : 28) ) , . . . 
bordl (:, i) ,[],[] , options) ; 


for j = 1 : length(xl ( : , 1) ) 

phi_modl ( j , i) = (pi/2) - atan(sqrt( ... 

(bmodl(2,i) + bmodl(8,i)*xl(j ,2) + ... 
bmodl(9,i)*xl(j ,3) + bmodl(10,i)*xl(j,4) + ... 
bmodl (11 , i)*xl (j ,5) + bmodl (12 , i) *xl ( j ,6) + ... 
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2*bmodl (23, i) *xl ( j , 1) ) ~2 + (bmodl(3,i) + ... 
bmodl (8, i)*xl (j , 1) + bmodl(13,i)*xl(j ,3) + ... 
bmodl(14,i)*xl(j ,4) + bmodl (15 , i) *xl ( j ,5) + ... 
bmodl (16 , i) *xl (j ,6) + 2*bmodl (24, i) *xl ( j ,2) ) ~2 + ... 
(bmodl (4, i) + bmodl (9 , i) *xl ( j , 1) + ... 
bmodl(13,i)*xl(j ,2) + bmodl (17, i) *xl ( j ,4) + ... 
bmodl(18,i)*xl(j ,5) + bmodl (19 , i) *xl ( j ,6) + ... 

2*bmodl (25, i) *xl ( j ,3) ) ~2 + (bmodl (5, i) + ... 
bmodl (10 , i) *xl (j , 1) + bmodl (14, i) *xl ( j ,2) + ... 
bmodl(17,i)*xl(j ,3) + bmodl (20, i)*xl (j ,5) + ... 
bmodl (21 , i)*xl (j ,6) + 2*bmodl (26 , i) *xl ( j ,4) ) ~2 + ... 
(bmodl (6, i) + bmodl (11 , i) *xl (j , 1) + ... 
bmodl (15, i)*xl (j ,2) + bmodl (18, i) *xl ( j ,3) + ... 
bmodl (20, i)*xl (j ,4) + bmodl (22, i)*xl (j ,6) + ... 

2*bmodl (27, i) *xl ( j , 5) ) ~2 + (bmodl (7, i) + ... 
bmodl (12, i)*xl (j , 1) + bmodl (16 , i) *xl ( j ,2) + ... 
bmodl (19, i)*xl (j ,3) + bmodl (21 , i)*xl (j ,4) + ... 
bmodl (22 , i) *xl (j ,5) + 2*bmodl (28, i) *xl ( j , 6) ) ~2) ) ; 

alpha_modl ( j , i) = (pi/2) + (1-gamma) *atan(sqrt ( ... 

(bmodl (2, i) + bmodl(8,i)*xl(j ,2) + ... 
bmodl(9,i)*xl(j ,3) + bmodl (10 , i) *xl (j ,4) + ... 
bmodl (11 , i)*xl (j ,5) + bmodl (12, i)*xl (j ,6) + ... 

2*bmodl (23, i) *xl ( j , 1) ) ~2 + (bmodl (3, i) + ... 
bmodl (8, i) *xl ( j , 1) + bmodl(13,i)*xl(j ,3) + ... 
bmodl(14,i)*xl(j ,4) + bmodl (15 , i) *xl ( j ,5) + ... 
bmodl (16, i)*xl (j ,6) + 2*bmodl (24, i) *xl ( j ,2) ) ~2 + ... 
(bmodl (4, i) + bmodl (9 , i) *xl ( j , 1) + ... 
bmodl(13,i)*xl(j ,2) + bmodl (17, i)*xl (j ,4) + ... 
bmodl(18,i)*xl(j ,5) + bmodl (19, i)*xl (j ,6) + ... 

2*bmodl (25, i) *xl ( j ,3) ) “2 + (bmodl (5, i) + ... 
bmodl (10, i)*xl (j , 1) + bmodl (14, i) *xl ( j ,2) + ... 
bmodl(17,i)*xl(j ,3) + bmodl (20, i)*xl (j ,5) + ... 
bmodl (21 , i)*xl (j ,6) + 2*bmodl (26 , i) *xl ( j ,4) ) ~2 + ... 
(bmodl (6, i) + bmodl (11 , i) *xl (j , 1) + ... 
bmodl (15, i)*xl (j ,2) + bmodl (18, i) *xl ( j ,3) + ... 
bmodl (20, i)*xl (j ,4) + bmodl (22 , i) *xl ( j ,6) + ... 

2*bmodl (27, i) *xl ( j , 5) ) ~2 + (bmodl (7, i) + ... 
bmodl (12, i)*xl (j , 1) + bmodl (16, i)*xl (j ,2) + ... 
bmodl (19, i)*xl (j ,3) + bmodl (21 , i) *xl ( j ,4) + ... 
bmodl (22 , i) *xl (j ,5) + 2*bmodl (28, i) *xl ( j , 6) ) ~2) ) ; 

ratio_modl ( j , i) = sin(phi_modl ( j , i) ) /sin(alpha_modl ( j , i) ) ; 
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end 

res_modl (1 : j , i) = (yobs_desl ( : , i) - xlm*bmodl ( : , i) ) ; 
for j = 1 : length(xl ( : , 1) ) 

dist_modl ( j , i) = ratio_raodl ( j , i) *res_modl ( j , i) ; 

end 

end 

for j = 1 : length(xl ( : , 1) ) 
dist_ordl ( j ,n+l) = 0; 

dist_raodl ( j ,n+l) = 0; 

dist_orthl ( j ,n+l) = 0; 

dist_ordl ( j ,n+2) = mean(dist_ordl (j , 1 :n) ) ; 

dist_modl ( j ,n+2) = mean(dist_modl (j , 1 :n) ) ; 

dist_orthl ( j ,n+2) = mean(dist_orthl ( j , 1 : n) ) ; 

dist_ordl ( j ,n+3) = var (dist_ordl ( j , 1 :n) ) ; 

dist_raodl ( j ,n+3) = var (dist_modl ( j , 1 :n) ) ; 

dist_orthl ( j ,n+3) = var (dist_orthl (j , 1 :n) ) ; 

dist_ordl ( j ,n+4) = dist_ordl ( j ,n+2) ~2 + dist_ordl (j ,n+3) ; 

dist_modl ( j ,n+4) = dist_modl ( j ,n+2) ~2 + dist_modl (j ,n+3) ; 

dist_orthl (j ,n+4) = dist_orthl (j ,n+2) "2 + dist_orthl (j ,n+3) ; 

end 

distl = [dist_ordl (1 : j ,n+2 :n+4) ; 000; ... 
dist_modl (1 : j ,n+2 :n+4) ; 000; ... 
dist_orthl(l : j ,n+2:n+4)] ; 


dist = [distl] ; 
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for i = l:n 

for k = 1 : N 

for j = 1 : length (xconf (:, 1) ) 
cal (i ,k) . conf ( j , : ) = ... 

mvnrnd( [xconf (j , 1) xconf (j, 2) xconf (j, 3) xconf (j, 4) ... 
xconf (j, 5) xconf (j, 6) yconf(j,l)], sigma); 
yobs (i) . conf (j ,k) = cal (i ,k) . conf (j ,7) ; 

cphi_ordl ( j , i) = (pi/2) - atan(sqrt( ... 

(bordl(2,i) + bordl (8 , i) *xconf (j , 2) + ... 
bordl (9 , i) *xconf ( j , 3) + bordl (10 , i) *xconf ( j ,4) + ... 
bordl (11 , i) *xconf ( j , 5) + bordl (12 ,i)*xconf (j ,6) + ... 
2*bordl (23 , i) *xconf (j , 1) ) ~2 + (bordl (3, i) + ... 
bordl (8, i) *xconf ( j , 1) + bordl (13 , i) *xconf (j , 3) + ... 
bordl(14,i)*xconf (j ,4) + bordl(15,i)*xconf (j ,5) + ... 
bordl (16 , i) *xconf ( j , 6) + ... 

2*bordl (24, i) *xconf (j , 2) ) ~2 + (bordl (4, i) + ... 
bordl (9 , i) *xconf ( j , 1) + bordl (13 , i) *xconf (j , 2) + ... 
bordl(17,i)*xconf (j ,4) + bordl (18 ,i)*xconf (j ,5) + ... 
bordl (19 , i) *xconf ( j , 6) + ... 

2*bordl (25 , i) *xconf (j , 3) ) ~2 + (bordl(5,i) + ... 
bordl (10 , i) *xconf (j , 1) + bordl (14, i) *xconf (j ,2) + ... 
bordl(17,i)*xconf (j ,3) + bordl (20, i)*xconf (j ,5) + ... 
bordl (21 , i) *xconf ( j , 6) + ... 

2*bordl (26 , i) *xconf (j ,4) ) ~2 + (bordl(6,i) + ... 
bordl (11 , i) *xconf (j , 1) + bordl(15,i)*xconf (j ,2) + ... 
bordl (18 , i) *xconf (j ,3) + bordl (20, i)*xconf (j ,4) + ... 
bordl (22 , i) *xconf ( j , 6) + ... 

2*bordl (27 , i) *xconf ( j , 5) ) ~2 + (bordl (7, i) + ... 
bordl (12 , i) *xconf (j , 1) + bordl(16,i)*xconf (j ,2) + ... 
bordl (19 , i) *xconf (j ,3) + bordl (21 , i) *xconf ( j ,4) + ... 
bordl (22 , i) *xconf ( j , 5) + 2*bordl (28 , i) *xconf ( j , 6) ) "2) ) ; 

cphi_orthl ( j , i) = (pi/2) - atan(sqrt( ... 

(borthl(2,i) + borthl (8, i) *xconf ( j ,2) + ... 
borthl (9 , i) *xconf (j ,3) + borthl (10 , i) *xconf ( j ,4) + ... 
borthl (11 , i) *xconf ( j ,5) + borthl (12 , i) *xconf ( j , 6) + ... 
2*borthl (23, i) *xconf ( j , 1) ) "2 + (borthl (3, i) + ... 
borthl (8 , i) *xconf (j , 1) + borthl (13, i) *xconf ( j , 3) + ... 
borthl (14, i) *xconf ( j ,4) + borthl (15 , i) *xconf ( j , 5) + ... 
borthl (16 , i) *xconf ( j , 6) + ... 

2*borthl(24,i)*xconf (j ,2)) “2 + (borthl (4, i) + ... 
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borthl (9, i)*xconf (j , 1) + borthl (13, i) *xconf ( j , 2) + ... 
borthl (17, i) *xconf ( j ,4) + borthl (18 , i) *xconf ( j , 5) + ... 
borthl (19 , i) *xconf ( j , 6) + ... 

2*borthl (25, i) *xconf ( j ,3) ) “2 + (borthl(5,i) + ... 
borthl (10 , i) *xconf (j , 1) + borthl (14, i)*xconf (j ,2) + ... 
borthl (17, i) *xconf ( j ,3) + borthl (20 , i) *xconf ( j , 5) + ... 
borthl (21 , i) *xconf ( j , 6) + ... 

2*borthl (26, i) *xconf ( j ,4) ) “2 + (borthl (6, i) + ... 
borthl (11 , i)*xconf (j , 1) + borthl (15, i)*xconf (j ,2) + ... 
borthl (18, i) *xconf ( j ,3) + borthl (20 , i) *xconf (j ,4) + ... 
borthl (22 , i) *xconf ( j , 6) + ... 

2*borthl (27, i) *xconf ( j , 5) ) “2 + (borthl (7, i) + ... 
borthl (12 , i) *xconf (j , 1) + borthl (16, i)*xconf (j ,2) + ... 
borthl (19, i)*xconf (j ,3) + borthl (21 , i) *xconf ( j ,4) + ... 
borthl (22 , i) *xconf ( j ,5) + ... 

2*borthl (28, i) *xconf ( j , 6) ) ~2) ) ; 

cphi_modl ( j , i) = (pi/2) - atan(sqrt( ... 

(bmodl(2,i) + bmodl (8 , i) *xconf (j , 2) + ... 
bmodl (9 , i) *xconf ( j , 3) + bmodl (10 , i) *xconf ( j ,4) + ... 
bmodl (11 , i) *xconf ( j , 5) + bmodl (12, i)*xconf (j ,6) + ... 
2*bmodl (23 , i) *xconf (j , 1) ) "2 + (bmodl (3, i) + ... 
bmodl (8, i) *xconf ( j , 1) + bmodl (13 , i) *xconf ( j , 3) + ... 
bmodl (14, i) *xconf (j ,4) + bmodl (15, i)*xconf (j ,5) + ... 
bmodl (16, i)*xconf (j ,6) + ... 

2*bmodl (24, i) *xconf (j , 2) ) ~2 + (bmodl (4, i) + ... 
bmodl (9, i)*xconf (j , 1) + bmodl (13 , i) *xconf (j , 2) + ... 
bmodl (17, i)*xconf (j ,4) + bmodl (18, i)*xconf (j ,5) + ... 
bmodl (19, i)*xconf (j ,6) + ... 

2*bmodl (25 , i) *xconf (j , 3) ) "2 + (bmodl (5, i) + ... 
bmodl (10, i)*xconf (j , 1) + bmodl (14, i) *xconf (j ,2) + ... 
bmodl (17, i)*xconf (j ,3) + bmodl (20, i)*xconf (j ,5) + ... 
bmodl (21 , i) *xconf ( j , 6) + ... 

2*bmodl (26 , i) *xconf (j ,4) ) "2 + (bmodl (6, i) + ... 
bmodl (11 , i) *xconf (j , 1) + bmodl (15, i)*xconf (j ,2) + ... 
bmodl (18, i)*xconf (j ,3) + bmodl (20, i)*xconf (j ,4) + ... 
bmodl (22, i)*xconf (j ,6) + ... 

2*bmodl (27 , i) *xconf (j , 5) ) ~2 + (bmodl (7, i) + ... 
bmodl (12, i)*xconf (j , 1) + bmodl (16, i)*xconf (j ,2) + ... 
bmodl(19,i)*xconf (j ,3) + bmodl (21 , i) *xconf ( j ,4) + ... 
bmodl (22, i)*xconf (j ,5) + ... 

2*bmodl (28 , i) *xconf ( j , 6) ) ~ 2 ) ) ; 

calpha_ordl (j , i) = (pi/2) - (gamma-1) * (atan(sqrt ( ... 
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(bordl(2,i) + bordl (8 , i) *xconf (j , 2) + ... 
bordl (9 , i) *xconf ( j , 3) + bordl (10 , i) *xconf ( j ,4) + ... 
bordl (11 , i) *xconf ( j , 5) + bordl(12,i)*xconf (j ,6) + ... 
2*bordl (23 , i) *xconf (j , 1) ) ~2 + (bordl(3,i) + ... 
bordl (8, i) *xconf ( j , 1) + bordl (13 , i) *xconf ( j , 3) + ... 
bordl (14, i) *xconf (j ,4) + bordl(15,i)*xconf (j ,5) + ... 
bordl (16, i)*xconf (j ,6) + ... 

2*bordl (24, i) *xconf (j , 2) ) ~2 + (bordl (4, i) + ... 
bordl (9 , i) *xconf ( j , 1) + bordl (13 , i) *xconf ( j , 2) + ... 
bordl(17,i)*xconf (j ,4) + bordl(18,i)*xconf (j ,5) + ... 
bordl (19, i)*xconf (j ,6) + ... 

2*bordl (25 , i) *xconf (j , 3) ) "2 + (bordl(5,i) + ... 
bordl (10, i)*xconf (j , 1) + bordl (14, i) *xconf (j ,2) + ... 
bordl (17, i)*xconf (j ,3) + bordl (20, i)*xconf (j ,5) + ... 
bordl (21 , i) *xconf ( j , 6) + ... 

2*bordl (26 , i) *xconf (j ,4) ) "2 + (bordl(6,i) + ... 
bordl (11 , i) *xconf (j , 1) + bordl(15,i)*xconf (j ,2) + ... 
bordl (18, i)*xconf (j ,3) + bordl (20, i)*xconf (j ,4) + ... 
bordl (22, i)*xconf (j ,6) + ... 

2*bordl(27,i)*xconf (j ,5))~2 + (bordl (7, i) + ... 
bordl (12, i)*xconf (j , 1) + bordl(16,i)*xconf (j ,2) + ... 
bordl (19, i)*xconf (j ,3) + bordl (21 , i) *xconf ( j ,4) + ... 
bordl (22, i)*xconf (j ,5) + ... 

2*bordl (28 , i) *xconf ( j , 6) ) ~2) ) ) ; 

calpha_modl (j , i) = (pi/2) - (gamma-1) * (atan(sqrt ( ... 
(bmodl(2,i) + bmodl (8 , i) *xconf (j , 2) + ... 
bmodl (9 , i) *xconf ( j , 3) + bmodl (10 , i) *xconf ( j ,4) + ... 
bmodl (11 , i)*xconf (j ,5) + bmodl (12, i)*xconf (j ,6) + ... 
2*bmodl (23 , i) *xconf (j , 1) ) "2 + (bmodl (3, i) + ... 
bmodl (8, i)*xconf (j , 1) + bmodl (13 , i) *xconf ( j , 3) + ... 
bmodl (14, i) *xconf (j ,4) + ... 

bmodl (15, i)*xconf (j ,5) + bmodl (16, i)*xconf (j ,6) + ... 
2*bmodl(24,i)*xconf (j ,2))~2 + (bmodl (4, i) + ... 
bmodl (9, i)*xconf (j , 1) + bmodl (13 , i) *xconf (j , 2) + ... 
bmodl (17, i)*xconf (j ,4) + bmodl (18, i)*xconf (j ,5) + ... 
bmodl (19, i)*xconf (j ,6) + ... 

2*bmodl (25 , i) *xconf (j , 3) ) ~2 + (bmodl (5, i) + ... 
bmodl (10, i)*xconf (j , 1) + bmodl (14, i) *xconf (j ,2) + ... 
bmodl (17, i)*xconf (j ,3) + bmodl (20, i)*xconf (j ,5) + ... 
bmodl (21 , i) *xconf ( j , 6) + ... 

2*bmodl (26 , i) *xconf (j ,4) ) ~2 + (bmodl (6, i) + ... 
bmodl (11 , i) *xconf (j , 1) + bmodl (15, i)*xconf (j ,2) + ... 
bmodl (18, i)*xconf (j ,3) + bmodl (20, i)*xconf (j ,4) + ... 
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bmodl (22, i)*xconf (j ,6) + ... 

2*bmodl (27 , i) *xconf (j , 5) ) "2 + (bmodl (7, i) + ... 
bmodl(12,i)*xconf (j , 1) + bmodl (16, i)*xconf (j ,2) + ... 
bmodl (19, i)*xconf (j ,3) + bmodl (21 , i) *xconf ( j ,4) + ... 
bmodl (22, i)*xconf (j ,5) + ... 

2*bmodl (28 , i) *xconf ( j , 6) ) ~2) ) ) ; 

calpha_orthl ( j , i) = (pi/2) - (gamma-1) * (atan(sqrt ( ... 
(borthl(2,i) + borthl (8, i) *xconf ( j ,2) + ... 
borthl(9,i)*xconf (j ,3) + borthl (10 , i) *xconf ( j ,4) + ... 
borthl (11 , i) *xconf ( j ,5) + borthl (12 , i) *xconf ( j , 6) + ... 
2*borthl (23, i) *xconf ( j , 1) ) ~2 + (borthl (3, i) + ... 
borthl (8, i)*xconf (j , 1) + borthl (13, i) *xconf ( j , 3) + ... 
borthl (14, i) *xconf ( j ,4) + borthl (15 , i) *xconf ( j , 5) + ... 
borthl (16 , i) *xconf ( j ,6) + ... 

2*borthl (24, i) *xconf ( j , 2) ) ~2 + (borthl (4, i) + ... 
borthl (9, i)*xconf (j , 1) + borthl (13, i) *xconf ( j , 2) + ... 
borthl (17, i) *xconf ( j ,4) + borthl (18 , i) *xconf ( j , 5) + ... 
borthl (19 , i) *xconf ( j , 6) + ... 

2*borthl (25, i) *xconf ( j ,3) ) ~2 + (borthl(5,i) + ... 
borthl (10 , i) *xconf (j , 1) + borthl(14,i)*xconf (j ,2) + ... 
borthl (17, i) *xconf ( j ,3) + borthl (20 , i) *xconf ( j , 5) + ... 
borthl (21 , i) *xconf ( j , 6) + ... 

2*borthl (26, i) *xconf ( j ,4) ) ~2 + (borthl(6,i) + ... 
borthl (11 , i) *xconf (j , 1) + borthl (15 , i) *xconf ( j , 2) + ... 
borthl (18, i) *xconf ( j ,3) + borthl (20 , i) *xconf (j ,4) + ... 
borthl (22 , i) *xconf ( j ,6) + ... 

2*borthl (27, i) *xconf ( j , 5) ) ~2 + (borthl (7, i) + ... 
borthl (12 , i) *xconf (j , 1) + borthl (16 , i) *xconf ( j , 2) + ... 
borthl (19 , i) *xconf ( j ,3) + borthl (21 , i) *xconf ( j ,4) + ... 
borthl (22 , i) *xconf ( j ,5) + ... 

2*borthl (28, i) *xconf ( j , 6) ) ~2) ) ) ; 

cratio_ordl (j ,k) = sin(cphi_ordl (j , i) )/ . . . 
sin(calpha_ordl ( j , i) ) ; 

cratio_modl (j ,k) = sin(cphi_modl (j , i) )/ . . . 
sin(calpha_modl ( j , i) ) ; 

cratio_orthl ( j ,k) = sin(cphi_orthl ( j , i) ) / . . . 
sin(calpha_orthl ( j , i) ) ; 


end 
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cres (i) . ordl (1 : j ,k) = (yobs (i) . conf ( : ,k) - ... 
xconfm*bordl ( : , i) ) ; 

cres (i) .modi (1 : j ,k) = (yobs (i) . conf (: ,k) - ... 
xconfm*bmodl ( : , i) ) ; 

cres (i) . orthl (1 : j ,k) = (yobs (i) . conf (: ,k) - ... 
xconfm*borthl ( : , i) ) ; 

for j = 1 : length (xconf (:, 1) ) 

cdist (i) . ordl ( j , i) = ... 

ratio_ordl ( j , i) *cres (i) . ordl ( j ,k) ; 

cdist (i) .modi (j , i) = ... 

ratio_modl ( j , i) *cres (i) .modi ( j ,k) ; 

cdist (i) . orthl (j , i) = ... 

ratio_orthl (j , i) *cres (i) . orthl (j ,k) ; 

end 

end 

for j = 1:38 

cdist (i) . ordl (j ,N+1) = 0; 
cdist(i) ,modl(j ,N+1) = 0; 
cdist (i) . orthl (j ,N+1) = 0; 

cdist (i) . ordl (j ,N+2) = mean(cdist (i) . ordl ( j , 1 : N) ) ; 
cdist (i) .modi (j ,N+2) = mean(cdist (i) .modi ( j , 1 : N) ) ; 
cdist(i) . orthl(j ,N+2) = mean(cdist (i) . orthl (j , 1 : N) ) ; 
cdist (i) . ordl (j ,N+3) = var (cdist (i) . ordl (j , 1 : N) ) ; 
cdist(i) ,modl(j ,N+3) = var (cdist (i) .modi (j , 1 : N) ) ; 
cdist(i) . orthl(j ,N+3) = var (cdist (i) . orthl (j , 1 : N) ) ; 


end 



141 


for j = 1:38 

ccdist_ordl ( j , i) = cdist (i) . ordl (j ,N+2) ; 
ccdist_modl ( j , i) = cdist (i) .modi (j ,N+2) ; 
ccdist_orthl ( j , i) = cdist (i) . orthl (j ,N+2) ; 

end 

end 

for j = 1:38 

ccdist_ordl (j ,n+l) = 0; 

ccdist_modl (j ,n+l) = 0; 

ccdist_orthl ( j ,n+l) = 0; 

ccdist_ordl (j ,n+2) = mean(ccdist_ordl ( j , 1 :n) ) ; 

ccdist_modl (j ,n+2) = mean(ccdist_modl ( j , 1 :n) ) ; 

ccdist_orthl ( j ,n+2) = mean(ccdist_orthl (j , 1 :n) ) ; 

ccdist_ordl (j ,n+3) = var (ccdist_ordl (j , 1 :n) ) ; 

ccdist_modl (j ,n+3) = var (ccdist_modl (j , 1 :n) ) ; 

ccdist_orthl ( j ,n+3) = var (ccdist_orthl (j , 1 :n) ) ; 

ccdist_ordl (j ,n+4) = ccdist_ordl ( j ,n+2) "2 + ccdist_ordl ( j ,n+3) ; 

ccdist_modl (j ,n+4) = ccdist_modl ( j ,n+2) "2 + ccdist_modl ( j ,n+3) ; 

ccdist_orthl ( j ,n+4) = ccdist_orthl ( j ,n+2) ~2 + ccdist_orthl ( j ,n+3) ; 

end 

ccdistl = [ccdist_ordl (1 : j ,n+2 :n+4) ; 000; ... 
ccdist_modl (1 : j ,n+2 : n+4) ; 000; ... 
ccdist_orthl (1 : j ,n+2 :n+4)] ; 

ccdist = [ccdistl] ; 


11 Mean and Variance of Coefficients over n Calibrations for 
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°/„°/o each Estimator and Design 

bmean = [mean(bordl (1 , : ) ) mean(bmodl (1 , : ) ) mean(borthl (1 , : ) ) ; ... 
mean(bordl (2 , : ) ) mean (bmodl (2 ,:) ) mean(borthl (2 , : ) ) ; ... 
mean(bordl (3, : ) ) mean(bmodl (3, : ) ) mean(borthl (3, : )) ; ... 
mean(bordl (4, : ) ) mean(bmodl (4, : ) ) mean(borthl (4, : )) ; ... 
mean(bordl (5 , : ) ) mean(bmodl (5 , : ) ) mean(borthl (5 , : ) ) ; ... 
mean(bordl (6 , : ) ) mean(bmodl (6 , : ) ) mean(borthl (6 , : ) ) ; ... 
mean(bordl (7, : ) ) mean(bmodl (7, : ) ) mean(borthl (7, : )) ; ... 
mean(bordl (8, : ) ) mean(bmodl (8, : ) ) mean(borthl (8, : )) ; ... 
mean(bordl (9 , : ) ) mean (bmodl (9 ,:) ) mean(borthl (9 , : ) ) ; ... 
mean(bordl (10 , : ) ) mean(bmodl (10 , : ) ) mean(borthl (10 , : ) ) ; ... 
mean(bordl (11 , : ) ) mean(bmodl (11 , : ) ) mean(borthl (11 , : ) ) ; ... 
mean(bordl (12 , : ) ) mean(bmodl (12 , : ) ) mean(borthl (12 , : ) ) ; ... 
mean(bordl (13 , : ) ) mean(bmodl (13, : ) ) mean(borthl (13, : ) ) ; ... 
mean(bordl (14, : ) ) mean(bmodl (14, : ) ) mean(borthl (14, : ) ) ; ... 
mean(bordl (15 , : ) ) mean(bmodl (15 , : ) ) mean(borthl (15, : ) ) ; ... 
mean(bordl (16 , : ) ) mean(bmodl (16 , : ) ) mean(borthl (16 , : ) ) ; ... 
mean(bordl (17 , : ) ) mean(bmodl (17, : ) ) mean(borthl (17, : ) ) ; ... 
mean(bordl (18 , : ) ) mean(bmodl (18, : ) ) mean(borthl (18, : ) ) ; ... 
mean(bordl (19 , : ) ) mean(bmodl (19 , : ) ) mean(borthl (19 , : ) ) ; ... 
mean(bordl (20 , : ) ) mean(bmodl (20 , : ) ) mean(borthl (20 , : ) ) ; ... 
mean(bordl (21 , : ) ) mean(bmodl (21 , : ) ) mean(borthl (21 , : ) ) ; ... 
mean(bordl (22 , : ) ) mean(bmodl (22 , : ) ) mean(borthl (22 , : ) ) ; ... 
mean(bordl (23 , : ) ) mean(bmodl (23, : ) ) mean(borthl (23, : ) ) ; ... 
mean(bordl (24, : ) ) mean(bmodl (24, : ) ) mean(borthl (24, : ) ) ; ... 
mean(bordl (25 , : ) ) mean(bmodl (25 , : ) ) mean(borthl (25 , : ) ) ; ... 
mean(bordl (26 , : ) ) mean(bmodl (26 , : ) ) mean(borthl (26 , : ) ) ; ... 
mean(bordl (27 , : ) ) mean(bmodl (27, : ) ) mean(borthl (27, : ) ) ; ... 
mean(bordl (28 , : ) ) mean(bmodl (28, : ) ) mean(borthl (28, : ) )] ; 

bvar = [var (bordl (1 , : ) ) var (bmodl (1 , : ) ) var (borthl (1 , : ) ) ; ... 
var(bordl(2, : )) var (bmodl (2, :) ) var (borthl (2 ,:)) ; ... 
var (bordl (3, : )) var (bmodl (3, :) ) var (borthl (3, :)) ; ... 
var (bordl (4, :) ) var (bmodl (4, :) ) var (borthl (4, :)) ; ... 
var (bordl (5, :) ) var (bmodl (5, :) ) var (borthl (5 ,:)) ; ... 
var (bordl (6, : )) var (bmodl (6, :) ) var (borthl (6 ,:)) ; ... 
var (bordl (7, :) ) var (bmodl (7, :) ) var (borthl (7, :)) ; ... 
var (bordl (8, : )) var (bmodl (8, :) ) var (borthl (8, :)) ; ... 
var (bordl (9, : )) var (bmodl (9, :) ) var (borthl (9 ,:)) ; ... 
var (bordl (10, : )) var (bmodl (10 ,:) ) var (borthl (10 ,:)) ; ... 
var (bordl (11 ,:) ) var (bmodl (11 ,:) ) var (borthl (11 ,:)) ; ... 
var (bordl (12, : )) var (bmodl (12 ,:) ) var (borthl (12 ,:)) ; ... 
var (bordl (13, :) ) var (bmodl (13, :) ) var (borthl (13, : )) ; ... 
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var (bordl (14, : ) ) 
var (bordl (15, : ) ) 
var(bordl(16, : )) 
var (bordl (17, : ) ) 
var (bordl (18, : ) ) 
var(bordl(19, : )) 
var (bordl (20, : )) 
var (bordl (21 , : ) ) 
var (bordl (22, : )) 
var (bordl (23, : ) ) 
var (bordl (24, : ) ) 
var (bordl (25, : )) 
var (bordl (26, : )) 
var (bordl (27, : ) ) 
var (bordl (28, : ) ) 


t = toe; 

tmin = t/60; 

disp( [’Elapsed time 


var (bmodl (14, : ) ) 
var (bmodl (15 , : ) ) 
var (bmodl (16 , : ) ) 
var (bmodl (17, : ) ) 
var (bmodl (18, : ) ) 
var (bmodl (19 , : ) ) 
var (bmodl (20 , : ) ) 
var (bmodl (21 , : ) ) 
var (bmodl (22 , : ) ) 
var (bmodl (23, : ) ) 
var (bmodl (24, : ) ) 
var (bmodl (25 , : ) ) 
var (bmodl (26 , : ) ) 
var (bmodl (27, : ) ) 
var (bmodl (28, : ) ) 


var (borthl ( 14 , : ) ) ; 
var (borthl ( 15 , : ) ) ; 
var (borthl ( 16 , : ) ) ; 
var (borthl ( 17 , : ) ) ; 
var (borthl ( 18 , : ) ) ; 
var (borthl ( 19 , : ) ) ; 
var (borthl (20 , : ) ) ; 
var (borthl (21 , : ) ) ; 
var (borthl (22 , : ) ) ; 
var (borthl (23, : )) ; 
var (borthl (24, : )) ; 
var (borthl (25, : )) ; 
var (borthl (26 , : ) ) ; 
var (borthl (27, : )) ; 
var (borthl (28, : ))] 


is ’, num2str (tmin) , ’ minutes.’]); 
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